5.1.1 Business Understanding

The BSPF includes the following steps of the the CRISP-DM Framework:

1. View the Business as a Machine
2. Understand the Drivers
3. Measure the Drivers
4. Uncover Problems & Opportunities

5.1.2 View The Business As A Machine

Review business_understanding.R file in 01_Business Understanding folder

Business Units: Department and Job Role Define Objectives: Retain high performers Assess Outcomes: TBD - Feedback from the business

Known - 16.1% attrition in the current state - is this a bad thing?

5.4.1 Group Features

Some features related to one and other.

It is always helpful to breakdown data collection activities in to strategic areas.

• Descriptive Features:
  • Age, DistanceFromHome, Gebder, MaritalStatus, NumComaniesWorked, Over18
• Employment Features:
  • Department, EmployeeCount, EmployeeNumber, JobInvolement, JobLevel, JobRole, JobSatistfaction
• Compensation Features:
  • DailyRate, HourlyRate, MonthlyIncome, MonthlyRate, PercentSalaryHike, StockOptionLevel
• Survey Results:
  • EnvironmentSatisfaction, JobSatisfaction, RelationshipSatisfaction, WorkLifeBalance
• Performance Data:
  • JobInvolvement, PerformanceRating
• Work-Life Balance
  • BusinessTravel, OverTime
• Training and Education
  • Education, EducationField, TrianingTimesLastYear
• Time-Based Features
  • TotalWorkingYears, YearsAtCompany, YeasInCurrentRole, YearsSinceLastPromotion, YearsWithCurrentManager

5.6.1 calculate_attrition_cost()

Create a function the mimics the Excel Calculator.

The function includes default values that equal the defaults that were introduced in Calculating the Cost of Turnover. So if we call the custom function calculate_attrition_cost with the default values, it returns 78483.33. If we change the default value of n = 1 (representing 1 employee) and change it to n = 200 then the function returns 15696666.67.

5.6.2 Calculating Cost By Job Level

Using code that we have already seen above and modified it slightly to focus on JobRole rather than Department, we can begin this part of the journey reviewing JobRole data. Note the we reuse a custom function - calculate_attrition_cost to return the cost_of_attrition column below.

dept_job_role_tbl %>% count(JobRole, Attrition) %>% group_by(JobRole) %>% mutate(pct = n/sum(n)) %>% ungroup() %>% 
     filter(Attrition %in% c("Yes")) %>% arrange(desc(pct)) %>% mutate(above_industry_avg = case_when(pct > 0.088 ~ "Yes", TRUE ~ "No")) %>% 
     mutate(cost_of_attrition = calculate_attrition_cost(n = n, salary = 80000))

## # A tibble: 9 x 6
##   JobRole        Attrition     n    pct above_industry_a~ cost_of_attriti~
##   <chr>          <chr>     <int>  <dbl> <chr>                        <dbl>
## 1 Sales Represe~ Yes          26 0.4    Yes                       2040567.
## 2 Human Resourc~ Yes          12 0.308  Yes                        941800.
## 3 Laboratory Te~ Yes          49 0.219  Yes                       3845683.
## 4 Sales Executi~ Yes          50 0.183  Yes                       3924167.
## 5 Research Scie~ Yes          43 0.166  Yes                       3374783.
## 6 Healthcare Re~ Yes           8 0.0762 No                         627867.
## 7 Manufacturing~ Yes           7 0.0569 No                         549383.
## 8 Manager        Yes           4 0.0449 No                         313933.
## 9 Research Dire~ Yes           2 0.0274 No                         156967.

Note the costs associated with Sales Executive and Laboratory Technician.

5.7.1 Streamlining The Counts

This is code commonly used:

dept_job_role_tbl %>% 
     
     group_by(JobRole, Attrition) %>%  summarise(n = n())

## # A tibble: 18 x 3
## # Groups:   JobRole [?]
##    JobRole                   Attrition     n
##    <chr>                     <chr>     <int>
##  1 Healthcare Representative No           97
##  2 Healthcare Representative Yes           8
##  3 Human Resources           No           27
##  4 Human Resources           Yes          12
##  5 Laboratory Technician     No          175
##  6 Laboratory Technician     Yes          49
##  7 Manager                   No           85
##  8 Manager                   Yes           4
##  9 Manufacturing Director    No          116
## 10 Manufacturing Director    Yes           7
## 11 Research Director         No           71
## 12 Research Director         Yes           2
## 13 Research Scientist        No          216
## 14 Research Scientist        Yes          43
## 15 Sales Executive           No          223
## 16 Sales Executive           Yes          50
## 17 Sales Representative      No           39
## 18 Sales Representative      Yes          26

The code above is replaced with this simpler code:

dept_job_role_tbl %>% count(JobRole, Attrition) 

## # A tibble: 18 x 3
##    JobRole                   Attrition     n
##    <chr>                     <chr>     <int>
##  1 Healthcare Representative No           97
##  2 Healthcare Representative Yes           8
##  3 Human Resources           No           27
##  4 Human Resources           Yes          12
##  5 Laboratory Technician     No          175
##  6 Laboratory Technician     Yes          49
##  7 Manager                   No           85
##  8 Manager                   Yes           4
##  9 Manufacturing Director    No          116
## 10 Manufacturing Director    Yes           7
## 11 Research Director         No           71
## 12 Research Director         Yes           2
## 13 Research Scientist        No          216
## 14 Research Scientist        Yes          43
## 15 Sales Executive           No          223
## 16 Sales Executive           Yes          50
## 17 Sales Representative      No           39
## 18 Sales Representative      Yes          26

count is part of dplyr!

Do not forget about the dplyr::tally function - very useful

tally() is a convenient wrapper for summarize that will either call n() or sum(n) depending on whether you’re tallying for the first time, or re-tallying. count() is similar but calls group_by() before and ungroup() after.

5.7.2 Streamlining The Count To Percentage Calculation

Pay attention - this is really good code using tidy evaluation!

We have this code so far in streamlining the attrition code:

dept_job_role_tbl %>% count(JobRole, Attrition) %>% group_by(JobRole) %>% mutate(pct = n/sum(n))

## # A tibble: 18 x 4
## # Groups:   JobRole [9]
##    JobRole                   Attrition     n    pct
##    <chr>                     <chr>     <int>  <dbl>
##  1 Healthcare Representative No           97 0.924 
##  2 Healthcare Representative Yes           8 0.0762
##  3 Human Resources           No           27 0.692 
##  4 Human Resources           Yes          12 0.308 
##  5 Laboratory Technician     No          175 0.781 
##  6 Laboratory Technician     Yes          49 0.219 
##  7 Manager                   No           85 0.955 
##  8 Manager                   Yes           4 0.0449
##  9 Manufacturing Director    No          116 0.943 
## 10 Manufacturing Director    Yes           7 0.0569
## 11 Research Director         No           71 0.973 
## 12 Research Director         Yes           2 0.0274
## 13 Research Scientist        No          216 0.834 
## 14 Research Scientist        Yes          43 0.166 
## 15 Sales Executive           No          223 0.817 
## 16 Sales Executive           Yes          50 0.183 
## 17 Sales Representative      No           39 0.6   
## 18 Sales Representative      Yes          26 0.4

Calculating percentages from counts is a common programming requirements. Why not create a function that generalizes this?

count_to_pct <- function(data, ..., col = n){
     
     grouping_vars_expr <- quos(...)
     col_expr <- enquo(col)
     
     ret <- data %>% group_by(!!! grouping_vars_expr) %>% 
          mutate(pct = (!! col_expr)/sum(!! col_expr)) %>% ungroup()
                 
     return(ret)     
}

• tidy functions always requires the data to take the first position of the function{#Tidy_Eval}
• ... enables passing multiple, un-named arguments to a function
  • Because the dots are not preselected, the user can flexibly add variables and the function will adapt!
• The 3rd function argument is a column name, in this case n since we are calculating a percentage from the counts in the n column

1. ... must be enclosed on quos or enquo for multiple columns or single columns, respectfully. The saves the column names as unevalutated expressions. They can be called later in the processing. Use quos for multiple columns and enquo for the single column we have called n
2. !!! called bang-bang-bang splices multiple grouping variables and evaluating them. (Use !! for single columns.)

Demo:

dept_job_role_tbl %>% count(Department,JobRole, Attrition) %>% count_to_pct(Department, JobRole)

## # A tibble: 21 x 5
##    Department             JobRole                   Attrition     n    pct
##    <chr>                  <chr>                     <chr>     <int>  <dbl>
##  1 Human Resources        Human Resources           No           27 0.692 
##  2 Human Resources        Human Resources           Yes          12 0.308 
##  3 Human Resources        Manager                   No           10 1     
##  4 Research & Development Healthcare Representative No           97 0.924 
##  5 Research & Development Healthcare Representative Yes           8 0.0762
##  6 Research & Development Laboratory Technician     No          175 0.781 
##  7 Research & Development Laboratory Technician     Yes          49 0.219 
##  8 Research & Development Manager                   No           46 0.958 
##  9 Research & Development Manager                   Yes           2 0.0417
## 10 Research & Development Manufacturing Director    No          116 0.943 
## # ... with 11 more rows

5.7.3 Streamlining The Attrition Assessment

Continuing the development of a streamlined function, review where we are so far:

dept_job_role_tbl %>% count(Department, JobRole, Attrition) %>% count_to_pct(Department, JobRole) %>% 
     filter(Attrition %in% c("Yes")) %>% arrange(desc(pct)) %>% 
     mutate(above_industry_avg = case_when(pct > 0.088 ~ "Yes", TRUE ~ "No"))

## # A tibble: 10 x 6
##    Department      JobRole         Attrition     n    pct above_industry_~
##    <chr>           <chr>           <chr>     <int>  <dbl> <chr>           
##  1 Sales           Sales Represen~ Yes          26 0.4    Yes             
##  2 Human Resources Human Resources Yes          12 0.308  Yes             
##  3 Research & Dev~ Laboratory Tec~ Yes          49 0.219  Yes             
##  4 Sales           Sales Executive Yes          50 0.183  Yes             
##  5 Research & Dev~ Research Scien~ Yes          43 0.166  Yes             
##  6 Research & Dev~ Healthcare Rep~ Yes           8 0.0762 No              
##  7 Sales           Manager         Yes           2 0.0645 No              
##  8 Research & Dev~ Manufacturing ~ Yes           7 0.0569 No              
##  9 Research & Dev~ Manager         Yes           2 0.0417 No              
## 10 Research & Dev~ Research Direc~ Yes           2 0.0274 No

Generalize filter - arrange - mutate statements:

assess_attrition <- function(data, attrition_col, attrition_value, baseline_pct){
     
     attrition_col_expr <- enquo(attrition_col)
     data %>% filter((!! attrition_col_expr) %in% attrition_value) %>% arrange(desc(pct)) %>% 
          mutate(above_industry_avg = case_when(pct > baseline_pct ~ "Yes", TRUE ~ "No"))}

• As always, tidy eval starts with a data (tibble)
• Next argument is the name of the Attrition column - it may not always be called Attrition
• 3rd argument are the values of the attrition_col - in our case the values are Yes and No.
• baseline_pct is hopefully self-explanatory!

1. Recall enquo saves a single column name as an expression that can be evaluated later. This way the name of the Attrition column can be whatever the user chooses!
2. pct does not change because that is created consistently in the count_to_pct function. (Of course this could be changed as an enhancement!)
3. Simple substitution for 0.088 to baseline_pct.

Use parenthesis to give Tidy Eval evaluations priority!

Function inputs in text format (attrition_value = “Yes”) do not require Tidy Eval

Demo the new code:

dept_job_role_tbl %>% 
     
     count(Department, JobRole, Attrition) %>% count_to_pct(Department, JobRole) %>% 
     assess_attrition(Attrition, attrition_value = "Yes", baseline_pct = 0.088) %>% #New code line
     mutate(cost_of_attrition = calculate_attrition_cost(n = n, salary = 80000))

## # A tibble: 10 x 7
##    Department JobRole Attrition     n    pct above_industry_~
##    <chr>      <chr>   <chr>     <int>  <dbl> <chr>           
##  1 Sales      Sales ~ Yes          26 0.4    Yes             
##  2 Human Res~ Human ~ Yes          12 0.308  Yes             
##  3 Research ~ Labora~ Yes          49 0.219  Yes             
##  4 Sales      Sales ~ Yes          50 0.183  Yes             
##  5 Research ~ Resear~ Yes          43 0.166  Yes             
##  6 Research ~ Health~ Yes           8 0.0762 No              
##  7 Sales      Manager Yes           2 0.0645 No              
##  8 Research ~ Manufa~ Yes           7 0.0569 No              
##  9 Research ~ Manager Yes           2 0.0417 No              
## 10 Research ~ Resear~ Yes           2 0.0274 No              
## # ... with 1 more variable: cost_of_attrition <dbl>

5.8.1 Data Manipulation For Visualization

Here is the code so far plus some data manipulation so we can produce the plot. {#fct_reorder}

dept_job_role_tbl %>% 
     
     count(Department, JobRole, Attrition) %>% count_to_pct(Department, JobRole) %>% 
     assess_attrition(Attrition, attrition_value = "Yes", baseline_pct = 0.088) %>% 
     mutate(cost_of_attrition = calculate_attrition_cost(n = n, salary = 80000)) %>% 

     # Data Manipulation 
     mutate(name = str_c(Department, JobRole, sep = ": ") %>%  as_factor()) %>% #LINE 1
     mutate(name = fct_reorder(name, cost_of_attrition)) %>% # LINE 2
     mutate(cost_text = str_c("$", format(cost_of_attrition / 1e6, digits = 2), "M", sep = "")) # LINE 3

## # A tibble: 10 x 9
##    Department JobRole Attrition     n    pct above_industry_~
##    <chr>      <chr>   <chr>     <int>  <dbl> <chr>           
##  1 Sales      Sales ~ Yes          26 0.4    Yes             
##  2 Human Res~ Human ~ Yes          12 0.308  Yes             
##  3 Research ~ Labora~ Yes          49 0.219  Yes             
##  4 Sales      Sales ~ Yes          50 0.183  Yes             
##  5 Research ~ Resear~ Yes          43 0.166  Yes             
##  6 Research ~ Health~ Yes           8 0.0762 No              
##  7 Sales      Manager Yes           2 0.0645 No              
##  8 Research ~ Manufa~ Yes           7 0.0569 No              
##  9 Research ~ Manager Yes           2 0.0417 No              
## 10 Research ~ Resear~ Yes           2 0.0274 No              
## # ... with 3 more variables: cost_of_attrition <dbl>, name <fct>,
## #   cost_text <chr>

• LINE 1: str_c from stringr concatenates Department and JobRole with a colon in between them.
  • factors are required in ggplot because they are used to order categorical variables (text will not work)
  • as_factor is part of forcats. This function creates levels in the order in which they appear, which will be the same on every platform. (Base R sorts in the current locale which can vary from place to place.)
• LINE 2: fct_reorder also from forcats reorders a factors numeric values by the magnitude of a different numeric value, in this case, cost_of_attrition.
• LINE 3: str_c concatenation again - pretty simple - but effective!

5.8.2 Visualization With ggplot2

OK, the data manipulation to plot the cost or attrition is complete. Now review the ggpolt code to produce the plot shown in the beginning on this section. It is easy - provided you do one step at a time!

dept_job_role_tbl %>% 
     
     count(Department, JobRole, Attrition) %>% count_to_pct(Department, JobRole) %>% 
     assess_attrition(Attrition, attrition_value = "Yes", baseline_pct = 0.088) %>% 
     mutate(cost_of_attrition = calculate_attrition_cost(n = n, salary = 80000)) %>% 

     # Data Manipulation 
     mutate(name = str_c(Department, JobRole, sep = ": ") %>%  as_factor()) %>% #must be factors to order
     mutate(name = fct_reorder(name, cost_of_attrition)) %>% # orders the names by cost
     mutate(cost_text = str_c("$", format(cost_of_attrition / 1e6, digits = 2), "M", sep = "")) %>% 
     
     # Plotting
     ggplot(aes(x = cost_of_attrition, y = name)) +
     geom_segment(aes(xend = 0, yend = name), color = palette_light()[[1]]) + # LINE 1
     geom_point(aes(size = cost_of_attrition), color = palette_light()[[1]]) + # LINE 2
     scale_x_continuous(label = scales::dollar) + # LINE 3
     geom_label(aes(label = cost_text, size = cost_of_attrition), hjust = "inward", color = palette_light()[[1]]) + # LINE 4
     theme_tq() + scale_size(name = c(4,5)) + # LINE 5
     labs(title = "Estimated Cost of Attrition: By Dept and Job Role", y = "", x = "Cost of Attrition") +
     theme(legend.position = "none")

• LINE 1: geom_segment adds lines to the plot canvas. The first ggplot statement locates the end of the line - at the cost_of_attrition. xend defines the start of the line segment at x = 0. yend make sure the line segment aligns with the correct y-axis label. Simple!
  • palette_light()[[1]] defines the color #2c3e5 from the tidyquant package.
  • the other colors in palette_light are: #2c3e50, #e31a1c, #18BC9C, #CCBE93, #a6cee3, #1f78b4, #b2df8a, #fb9a99, #fdbf6f, #ff7f00, #cab2d6, #6a3d9a
• LINE 2: size = cost_of_attrition adjusts the size of the segment endpoint to the attrition cost magnitude (could have included this in the first aes statement)
• LINE 3: Formatting the x-labels with the scales package to get the values in dollar format. Remember the scale_x_continuous ggplot function!
• LINE 4: geom_label simply adds labels to points on the plot using cost_text` that was created easier in the plot data manipulation.
  • Cool the label size can be dynamically set by the attrition cost.
• LINE 5: theme_tq from tidyquant is a general theme that adjusts background, faceting palette, legend position and more. Works well for conservative business publications.
  • scale_size adjusts the max and min size of elements to prevent large/small values from becoming too large/small. The label sizes are now limited to size = 3 making the lower values label text now readable. Cool!

6.1.1 Character Data

skimr does provide a good overview of the data automatically separating character and numeric data. For character data, it is simple to see if any features have more than perhaps 9 unique features, it might suggest binning. The numeric data with histograms makes it easy to see which variable might be categorical data very quickly.

skim(train_raw_tbl)

## Skim summary statistics
##  n obs: 1250 
##  n variables: 35 
## 
## -- Variable type:character ---------------------------------------------------------------------------------------------
##        variable missing complete    n min max empty n_unique
##       Attrition       0     1250 1250   2   3     0        2
##  BusinessTravel       0     1250 1250  10  17     0        3
##      Department       0     1250 1250   5  22     0        3
##  EducationField       0     1250 1250   5  16     0        6
##          Gender       0     1250 1250   4   6     0        2
##         JobRole       0     1250 1250   7  25     0        9
##   MaritalStatus       0     1250 1250   6   8     0        3
##          Over18       0     1250 1250   1   1     0        1
##        OverTime       0     1250 1250   2   3     0        2
## 
## -- Variable type:numeric -----------------------------------------------------------------------------------------------
##                  variable missing complete    n     mean      sd   p0
##                       Age       0     1250 1250    36.95    9.08   18
##                 DailyRate       0     1250 1250   799.74  403.52  102
##          DistanceFromHome       0     1250 1250     9.1     8.1     1
##                 Education       0     1250 1250     2.91    1.03    1
##             EmployeeCount       0     1250 1250     1       0       1
##            EmployeeNumber       0     1250 1250  1022.71  600.42    1
##   EnvironmentSatisfaction       0     1250 1250     2.71    1.09    1
##                HourlyRate       0     1250 1250    65.82   20.14   30
##            JobInvolvement       0     1250 1250     2.75    0.71    1
##                  JobLevel       0     1250 1250     2.08    1.12    1
##           JobSatisfaction       0     1250 1250     2.73    1.11    1
##             MonthlyIncome       0     1250 1250  6569.4  4777.89 1009
##               MonthlyRate       0     1250 1250 14294.5  7071.61 2097
##        NumCompaniesWorked       0     1250 1250     2.71    2.49    0
##         PercentSalaryHike       0     1250 1250    15.23    3.69   11
##         PerformanceRating       0     1250 1250     3.16    0.36    3
##  RelationshipSatisfaction       0     1250 1250     2.71    1.09    1
##             StandardHours       0     1250 1250    80       0      80
##          StockOptionLevel       0     1250 1250     0.8     0.86    0
##         TotalWorkingYears       0     1250 1250    11.4     7.79    0
##     TrainingTimesLastYear       0     1250 1250     2.81    1.29    0
##           WorkLifeBalance       0     1250 1250     2.76    0.71    1
##            YearsAtCompany       0     1250 1250     7.09    6.21    0
##        YearsInCurrentRole       0     1250 1250     4.26    3.66    0
##   YearsSinceLastPromotion       0     1250 1250     2.21    3.2     0
##      YearsWithCurrManager       0     1250 1250     4.16    3.59    0
##      p25     p50      p75  p100     hist
##    30       36      43       60 <U+2582><U+2583><U+2586><U+2587><U+2585><U+2583><U+2582><U+2582>
##   464      798    1156     1498 <U+2587><U+2587><U+2587><U+2587><U+2587><U+2587><U+2587><U+2587>
##     2        7      13.75    29 <U+2587><U+2585><U+2583><U+2581><U+2581><U+2581><U+2581><U+2581>
##     2        3       4        5 <U+2582><U+2585><U+2581><U+2587><U+2581><U+2586><U+2581><U+2581>
##     1        1       1        1 <U+2581><U+2581><U+2581><U+2587><U+2581><U+2581><U+2581><U+2581>
##   486.25  1020.5  1553.75  2068 <U+2587><U+2587><U+2587><U+2587><U+2587><U+2587><U+2587><U+2587>
##     2        3       4        4 <U+2585><U+2581><U+2585><U+2581><U+2581><U+2587><U+2581><U+2587>
##    48       66      83      100 <U+2586><U+2587><U+2587><U+2587><U+2586><U+2587><U+2587><U+2587>
##     2        3       3        4 <U+2581><U+2581><U+2583><U+2581><U+2581><U+2587><U+2581><U+2582>
##     1        2       3        5 <U+2587><U+2587><U+2581><U+2583><U+2581><U+2582><U+2581><U+2581>
##     2        3       4        4 <U+2585><U+2581><U+2585><U+2581><U+2581><U+2587><U+2581><U+2587>
##  2935.25  4903.5  8395    19999 <U+2587><U+2587><U+2583><U+2582><U+2581><U+2581><U+2581><U+2582>
##  8191.25 14328   20337.25 26999 <U+2587><U+2587><U+2587><U+2587><U+2587><U+2587><U+2587><U+2586>
##     1        2       4        9 <U+2587><U+2582><U+2582><U+2582><U+2581><U+2581><U+2581><U+2581>
##    12       14      18       25 <U+2587><U+2587><U+2583><U+2583><U+2582><U+2582><U+2582><U+2581>
##     3        3       3        4 <U+2587><U+2581><U+2581><U+2581><U+2581><U+2581><U+2581><U+2582>
##     2        3       4        4 <U+2585><U+2581><U+2586><U+2581><U+2581><U+2587><U+2581><U+2587>
##    80       80      80       80 <U+2581><U+2581><U+2581><U+2587><U+2581><U+2581><U+2581><U+2581>
##     0        1       1        3 <U+2587><U+2581><U+2587><U+2581><U+2581><U+2582><U+2581><U+2581>
##     6       10      15       40 <U+2583><U+2587><U+2582><U+2582><U+2582><U+2581><U+2581><U+2581>
##     2        3       3        6 <U+2581><U+2581><U+2587><U+2587><U+2581><U+2582><U+2582><U+2581>
##     2        3       3        4 <U+2581><U+2581><U+2583><U+2581><U+2581><U+2587><U+2581><U+2582>
##     3        5      10       40 <U+2587><U+2585><U+2581><U+2581><U+2581><U+2581><U+2581><U+2581>
##     2        3       7       18 <U+2587><U+2583><U+2581><U+2585><U+2581><U+2581><U+2581><U+2581>
##     0        1       3       15 <U+2587><U+2582><U+2581><U+2581><U+2581><U+2581><U+2581><U+2581>
##     2        3       7       17 <U+2587><U+2583><U+2581><U+2583><U+2581><U+2581><U+2581><U+2581>

6.1.2 Character Data

train_raw_tbl %>% select_if(is.character) %>% glimpse()

## Observations: 1,250
## Variables: 9
## $ Attrition      <chr> "Yes", "No", "Yes", "No", "No", "No", "No", "No...
## $ BusinessTravel <chr> "Travel_Rarely", "Travel_Frequently", "Travel_R...
## $ Department     <chr> "Sales", "Research & Development", "Research & ...
## $ EducationField <chr> "Life Sciences", "Life Sciences", "Other", "Lif...
## $ Gender         <chr> "Female", "Male", "Male", "Female", "Male", "Ma...
## $ JobRole        <chr> "Sales Executive", "Research Scientist", "Labor...
## $ MaritalStatus  <chr> "Single", "Married", "Single", "Married", "Marr...
## $ Over18         <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y...
## $ OverTime       <chr> "Yes", "No", "Yes", "Yes", "No", "No", "Yes", "...

Character data does not have levels so use unique instead. map from purrr maps unique across each feature.

train_raw_tbl %>% select_if(is.character) %>% map(unique) %>% head()

## $Attrition
## [1] "Yes" "No" 
## 
## $BusinessTravel
## [1] "Travel_Rarely"     "Travel_Frequently" "Non-Travel"       
## 
## $Department
## [1] "Sales"                  "Research & Development"
## [3] "Human Resources"       
## 
## $EducationField
## [1] "Life Sciences"    "Other"            "Medical"         
## [4] "Marketing"        "Technical Degree" "Human Resources" 
## 
## $Gender
## [1] "Female" "Male"  
## 
## $JobRole
## [1] "Sales Executive"           "Research Scientist"       
## [3] "Laboratory Technician"     "Manufacturing Director"   
## [5] "Healthcare Representative" "Manager"                  
## [7] "Sales Representative"      "Research Director"        
## [9] "Human Resources"

This improves upon unique by including the counts for each character feature.

train_raw_tbl %>% select_if(is.character) %>% map(table) %>% head()

## $Attrition
## 
##   No  Yes 
## 1049  201 
## 
## $BusinessTravel
## 
##        Non-Travel Travel_Frequently     Travel_Rarely 
##               124               239               887 
## 
## $Department
## 
##        Human Resources Research & Development                  Sales 
##                     49                    832                    369 
## 
## $EducationField
## 
##  Human Resources    Life Sciences        Marketing          Medical 
##               23              520              135              392 
##            Other Technical Degree 
##               66              114 
## 
## $Gender
## 
## Female   Male 
##    491    759 
## 
## $JobRole
## 
## Healthcare Representative           Human Resources 
##                       105                        39 
##     Laboratory Technician                   Manager 
##                       224                        89 
##    Manufacturing Director         Research Director 
##                       123                        73 
##        Research Scientist           Sales Executive 
##                       259                       273 
##      Sales Representative 
##                        65

Counts are nice but proportions are better! Below, anonymous function is used. {#anonymous}

The ~ indicates a formula starts

Also understand that train_raw_tbl %>% select_if(is.character) %>% map(~ table(.)) is the sames as train_raw_tbl %>% select_if(is.character) %>% map(table).

# Anonymous function with ~  
train_raw_tbl %>% select_if(is.character) %>% map(~ table(.) %>% prop.table()) 

## $Attrition
## .
##   No  Yes 
## 0.84 0.16 
## 
## $BusinessTravel
## .
##        Non-Travel Travel_Frequently     Travel_Rarely 
##             0.099             0.191             0.710 
## 
## $Department
## .
##        Human Resources Research & Development                  Sales 
##                  0.039                  0.666                  0.295 
## 
## $EducationField
## .
##  Human Resources    Life Sciences        Marketing          Medical 
##            0.018            0.416            0.108            0.314 
##            Other Technical Degree 
##            0.053            0.091 
## 
## $Gender
## .
## Female   Male 
##   0.39   0.61 
## 
## $JobRole
## .
## Healthcare Representative           Human Resources 
##                     0.084                     0.031 
##     Laboratory Technician                   Manager 
##                     0.179                     0.071 
##    Manufacturing Director         Research Director 
##                     0.098                     0.058 
##        Research Scientist           Sales Executive 
##                     0.207                     0.218 
##      Sales Representative 
##                     0.052 
## 
## $MaritalStatus
## .
## Divorced  Married   Single 
##     0.22     0.46     0.32 
## 
## $Over18
## .
## Y 
## 1 
## 
## $OverTime
## .
##   No  Yes 
## 0.72 0.28

6.1.3 Numeric Data

map_df provides a much cleaner output! Outputs a dataframe rather than a list.

train_raw_tbl %>% select_if(is.numeric) %>% map_df(~ unique(.) %>% length())

## # A tibble: 1 x 26
##     Age DailyRate DistanceFromHome Education EmployeeCount EmployeeNumber
##   <int>     <int>            <int>     <int>         <int>          <int>
## 1    43       808               29         5             1           1250
## # ... with 20 more variables: EnvironmentSatisfaction <int>,
## #   HourlyRate <int>, JobInvolvement <int>, JobLevel <int>,
## #   JobSatisfaction <int>, MonthlyIncome <int>, MonthlyRate <int>,
## #   NumCompaniesWorked <int>, PercentSalaryHike <int>,
## #   PerformanceRating <int>, RelationshipSatisfaction <int>,
## #   StandardHours <int>, StockOptionLevel <int>, TotalWorkingYears <int>,
## #   TrainingTimesLastYear <int>, WorkLifeBalance <int>,
## #   YearsAtCompany <int>, YearsInCurrentRole <int>,
## #   YearsSinceLastPromotion <int>, YearsWithCurrManager <int>

{#gather_example}

train_raw_tbl %>% select_if(is.numeric) %>% map_df(~ unique(.) %>% length()) %>% 
     gather() %>% arrange(value) %>% filter(value <= 10) # features with more than 10 might be nonessential or categorical

## # A tibble: 13 x 2
##    key                      value
##    <chr>                    <int>
##  1 EmployeeCount                1
##  2 StandardHours                1
##  3 PerformanceRating            2
##  4 EnvironmentSatisfaction      4
##  5 JobInvolvement               4
##  6 JobSatisfaction              4
##  7 RelationshipSatisfaction     4
##  8 StockOptionLevel             4
##  9 WorkLifeBalance              4
## 10 Education                    5
## 11 JobLevel                     5
## 12 TrainingTimesLastYear        7
## 13 NumCompaniesWorked          10

Numeric features that have fewer levels are likely to be discrete and numeric variables with more are likely continuous. (Note above the features where value = 1 may indicate non-essential or zero-variance variables.) Above are features that might be converted to factors. The table below suggests the features remain continuous.

train_raw_tbl %>% select_if(is.numeric) %>% map_df(~ unique(.) %>% length()) %>% 
     gather() %>% arrange(value) %>% filter(value > 10) # features with more than 10 might be nonessential or categorical

## # A tibble: 13 x 2
##    key                     value
##    <chr>                   <int>
##  1 PercentSalaryHike          15
##  2 YearsSinceLastPromotion    16
##  3 YearsWithCurrManager       18
##  4 YearsInCurrentRole         19
##  5 DistanceFromHome           29
##  6 YearsAtCompany             37
##  7 TotalWorkingYears          40
##  8 Age                        43
##  9 HourlyRate                 71
## 10 DailyRate                 808
## 11 MonthlyIncome            1161
## 12 MonthlyRate              1217
## 13 EmployeeNumber           1250

6.1.4 Visualizing The Feature-Target Interactions

To begin the visualization of the data, select the descriptive variables (serves as an example). There is a mix of character and continuous features.

train_raw_tbl %>% select(Attrition, Age, Gender, MaritalStatus, NumCompaniesWorked, Over18, DistanceFromHome)

## # A tibble: 1,250 x 7
##    Attrition   Age Gender MaritalStatus NumCompaniesWor~ Over18
##    <chr>     <dbl> <chr>  <chr>                    <dbl> <chr> 
##  1 Yes          41 Female Single                       8 Y     
##  2 No           49 Male   Married                      1 Y     
##  3 Yes          37 Male   Single                       6 Y     
##  4 No           33 Female Married                      1 Y     
##  5 No           27 Male   Married                      9 Y     
##  6 No           32 Male   Single                       0 Y     
##  7 No           59 Female Married                      4 Y     
##  8 No           30 Male   Divorced                     1 Y     
##  9 No           38 Male   Single                       0 Y     
## 10 No           36 Male   Married                      6 Y     
## # ... with 1,240 more rows, and 1 more variable: DistanceFromHome <dbl>

Apply ggpair function:

train_raw_tbl %>% select(Attrition, Age, Gender, MaritalStatus, NumCompaniesWorked, Over18, DistanceFromHome) %>% ggpairs()

Customize the output to make it a bit easier to interpret.

• It colors the Attrition to make it stand out
• Removes the lower triangle - it is not quite as helpful or understandable
• customizes what is shown on the diagonal and changes alpha

train_raw_tbl %>% select(Attrition, Age, Gender, MaritalStatus, NumCompaniesWorked, Over18, DistanceFromHome) %>% 
     ggpairs(aes(color = Attrition), lower = "blank", legend = 1, diag = list(continuous = wrap("densityDiag", alpha = 0.5))) +
     theme(legend.position = "bottom")

#ggsave("./images/plot_ggpairs.jpg")

6.1.5 Custom Function: plot_ggpairs

Take the visualization above and create a function.

• color= NULL will output normal grey. However, best to select a feature to highlight (typically the label)
  • rlang::quo_is_null simply tests for if color_expr == NULL. Test by simply typing color_expr in the console
• Because color_name is a string after using quo_name, aes_string must be used

plot_ggpairs <- function(data, color = NULL, density_alpha = 0.5){
     color_expr <- enquo(color)
     if(rlang::quo_is_null(color_expr)){gdata %>% ggpairs(lowe = "blank")
     } else {
          color_name <- quo_name(color_expr)
          g <- data %>%
               ggpairs(mapping = aes_string(color = color_name),
                       lower = "blank", legend = 1, diag = list(continuous = wrap("densityDiag", 
                                   alpha = density_alpha))) + theme(legend.position = "bottom")}
     return(g)}

6.1.6 Evaluate ggpairs Output

Evaluate the data features by category and understand the importance of the features using the new ggpairs function. How great is that?

1. Descriptive Features: age, gender, marital status, etc

train_raw_tbl  %>% select(Attrition, Age, Gender, MaritalStatus, NumCompaniesWorked, Over18, DistanceFromHome) %>% plot_ggpairs(Attrition)

Observations

• AgeYounger people seem to leave more than older workers.
  
• DistanceFromHome appears skewed - further away leave more frequently
• Num Companies Worked appears more skewed for workers that have had more jobs
• Marital Status interesting: No goes small large medium where Yes goes small medium large (height of histogram bars)

2. Employment Features

train_raw_tbl %>% select(Attrition, contains("employee"), contains("department"), contains("job")) %>% plot_ggpairs(Attrition)

Contains is NOT case sensitive. Also provides end_with option

Observations

• EmployeeCount and EmployeeNumber are not useful - will not be used in model
• Department: Higher proportion evident in Dept-Dept histogram in far right histogram bar
  • Recall there only only 3 departments in this data (HR, Research and Sales)
• JobInvolvement: Density Plot shows higher 3rd spike and 4th to less extent - more likely to stay
• JobLevel: Density plot shows deviation
• JobRoles - some roles appear to have higher attrition rate by evaluating the histogram
• JobStatisfaction: 1 - higher density of attrition than others

3. Compensation Features

train_raw_tbl %>% select(Attrition, contains("income"), contains("rate"), contains("salary"), 
                         contains("stock")) %>% plot_ggpairs(Attrition)

Observations

• MonthlyRate: Those that are leaving have a lower Monthly Income
• PercentSalaryHike: It’s difficult to deduce anything based on the visualization
• StockOptionLevel: Those that are staying have a higher stock option level

4. Survey Results

train_raw_tbl %>% select(Attrition, contains("satisfaction"), contains("life")) %>% plot_ggpairs(Attrition)

Observations

• EnvironmentSatisfaction: A higher proportion of those leaving have a low environment satisfaction level
• WorkLifeBalance: Those that are staying have a higher density of 2’s and 3’s

5. Performance Data

train_raw_tbl %>% select(Attrition, contains("performance"), contains("involvement")) %>% plot_ggpairs(Attrition)

• Observations:
• JobInvolvement: Those that are leaving have a lower density of 3’s and 4’s

6. Work-Life Features

train_raw_tbl %>% select(Attrition, contains("overtime"), contains("travel")) %>% plot_ggpairs(Attrition)

Observations

• Overtime: The proportion of those leaving that are working Over Time are high compared to those that are not leaving

7. Training and Education

train_raw_tbl %>% select(Attrition, contains("training"), contains("education")) %>% plot_ggpairs(Attrition)

Observations

• TrainingTimesLastYear: People that leave tend to have less annual training

8. Time-Based Features

train_raw_tbl %>% select(Attrition, contains("years")) %>% plot_ggpairs(Attrition)

Observations

• YearsAtCompany: People that leave tend to have less working years at the company
• YearsSinceLastPromotion: It’s difficult to deduce anything based on the visualization

7.1.1 Tidying The Data

Explore the data definitions file - note that many of the features are coded - numbers used as values for descriptions for categorical data.

path_data_definitions <- "00_Data/telco_data_definitions.xlsx"
definitions_raw_tbl <- read_excel(path_data_definitions, sheet = 1, col_names = FALSE)
definitions_raw_tbl

## # A tibble: 35 x 2
##    X__1                    X__2             
##    <chr>                   <chr>            
##  1 Education               1 'Below College'
##  2 <NA>                    2 'College'      
##  3 <NA>                    3 'Bachelor'     
##  4 <NA>                    4 'Master'       
##  5 <NA>                    5 'Doctor'       
##  6 <NA>                    <NA>             
##  7 EnvironmentSatisfaction 1 'Low'          
##  8 <NA>                    2 'Medium'       
##  9 <NA>                    3 'High'         
## 10 <NA>                    4 'Very High'    
## # ... with 25 more rows

X__1 = Feature Name; x__2 = Feature Code + Feature Description

Tidying this data will allow us to replace numbers with something more meaningful using tidyverse functions.

• fill replaces missing values (NAs) with the closet entry (previous if .direction = "down" or next if .direction = "up")
  • For example, this will fill NAs in Education under X__1 with Education - very cool!
• Remove NAs under X__2
• Break X__2 into 2 columns: one with the number, the other with the description
  • separate is from tidyr
  • remove simply removes the original X__2 column
  • Need to start using this!
• Rename columns to make them more friendly - X__1 does not mean anything!
• Remove the trailing tick mark from the new value column
  • str_replace replaces a matched pattern in a string with a replacement argument

 definitions_tbl <- definitions_raw_tbl %>% fill(X__1, .direction = "down") %>% filter(!is.na(X__2)) %>% 
     separate(X__2, into = c("key", "value"), sep = " '", remove = TRUE) %>% 
     rename(column_name = X__1) %>% mutate(key = as.numeric(key)) %>%  
     mutate(value = value %>% str_replace(pattern = "'", replacement = ""))
definitions_tbl

## # A tibble: 29 x 3
##    column_name               key value        
##    <chr>                   <dbl> <chr>        
##  1 Education                   1 Below College
##  2 Education                   2 College      
##  3 Education                   3 Bachelor     
##  4 Education                   4 Master       
##  5 Education                   5 Doctor       
##  6 EnvironmentSatisfaction     1 Low          
##  7 EnvironmentSatisfaction     2 Medium       
##  8 EnvironmentSatisfaction     3 High         
##  9 EnvironmentSatisfaction     4 Very High    
## 10 JobInvolvement              1 Low          
## # ... with 19 more rows

7.1.2 Mapping Over Lists

definitions_tbl has multiple datasets in one dataframe. They will need to be integrated separately into the training and testing datasets. This will make more sense in a bit. This process starts by making definitions_tbl into a list.

This concept is new - read and study carefully!

• split splits a dataframe into multiple dataframes contained within a list. Supply a column name as a vector (meaning .$column_name). Therefore, Education, EnvirnmentSatisfaction and JobInvolvement will be separated into individual dataframes. Example:

definitions_list[[3]]
# A tibble: 4 x 3
  column_name      key value    
  <chr>          <dbl> <chr>    
1 JobInvolvement     1 Low      
2 JobInvolvement     2 Medium   
3 JobInvolvement     3 High     
4 JobInvolvement     4 Very High

• Use map and an anonymous function (the ~) to remove `column_name from each of the dataframes stored in the newly created list
  • Example:

  > definitions_list[[3]]
  # A tibble: 4 x 2
    key value    
  <dbl> <chr>    
  1     1 Low      
  2     2 Medium   
  3     3 High     
  4     4 Very High

• Use map again with an anonymous function to each entry in the column value into a factor

definitions_list <- definitions_tbl %>% 
     split(.$column_name) %>% #split df into multiple dfs within a list
     map(~ select(., -column_name)) %>% # remove column_name col leaving key and value
     map(~ mutate(., value = as_factor(value)))# Creates factors in the order they appear - not alphabetically like as.factor
definitions_list[[3]]

## # A tibble: 4 x 2
##     key value    
##   <dbl> <fct>    
## 1     1 Low      
## 2     2 Medium   
## 3     3 High     
## 4     4 Very High

IMPORTANT: forcats::as_factor is different from base::as.factor. The former creates factor in the order in which they appear as opposed to alphabetically with base::as.factor

• Rename the names from key and value to something like Education and Education_value
  • first get the name of the dataframe in the list and assign it the variable list_name
  • Then set the column names using list_name and making the value column into list_name pasted with _value - really not that hard!

for(i in seq_along(definitions_list)){
     list_name <- names(definitions_list)[i]
     colnames(definitions_list[[i]]) <- c(list_name, paste0(list_name, "_value"))}
definitions_list[[3]]

## # A tibble: 4 x 2
##   JobInvolvement JobInvolvement_value
##            <dbl> <fct>               
## 1              1 Low                 
## 2              2 Medium              
## 3              3 High                
## 4              4 Very High

This can now be merged with the training dataset!

Iteratively join the dataframes within the definitions list with the main dataframe (training data)

This is new - study and learn!

• Make the train_raw_tbl at list and call it HRData - weird, huh!
  • make a list but the original training DF is the only entry called HRData
• append adds elements to a vector or list
  • The list grows because the first element is the training DF followed by all the elements in the definitions list
  • everything is in one list - helpful to use lists to collect objects that need to be iterated over with purrr
  • recall we are doing this so values are replaced by descriptions in the training data
• Use reduce which iteratively applies a user specified function to successive binary sets of objects. A 3-element vector would have a function applied to the first 2 elements and that output would then have the function applied with the 3rd element
  • here it adds columns like EducationField joined on the Education numerical value - human readable!
• Remove the original training data field names, like Education since we have just joined the more descriptive fields like Education_value
  • Using negative one_of to remove them
• Rename the remaining fields like Education_value to a more useful nomenclature
  • using set_names from base to do this leverage str_replace_all to remove _value
• Resort the dataframe

data_merged_tbl <- list(HRData = train_raw_tbl) %>% append(definitions_list, after = 1) %>% 
     reduce(left_join) %>% 
     select(-one_of(names(definitions_list))) %>% 
     set_names(str_replace_all(names(.), pattern = "_value", replacement = "")) %>% 
     select(sort(names(.)))
glimpse(data_merged_tbl)

## Observations: 1,250
## Variables: 35
## $ Age                      <dbl> 41, 49, 37, 33, 27, 32, 59, 30, 38, 3...
## $ Attrition                <chr> "Yes", "No", "Yes", "No", "No", "No",...
## $ BusinessTravel           <chr> "Travel_Rarely", "Travel_Frequently",...
## $ DailyRate                <dbl> 1102, 279, 1373, 1392, 591, 1005, 132...
## $ Department               <chr> "Sales", "Research & Development", "R...
## $ DistanceFromHome         <dbl> 1, 8, 2, 3, 2, 2, 3, 24, 23, 27, 16, ...
## $ Education                <fct> College, Below College, College, Mast...
## $ EducationField           <chr> "Life Sciences", "Life Sciences", "Ot...
## $ EmployeeCount            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ EmployeeNumber           <dbl> 1, 2, 4, 5, 7, 8, 10, 11, 12, 13, 14,...
## $ EnvironmentSatisfaction  <fct> Medium, High, Very High, Very High, L...
## $ Gender                   <chr> "Female", "Male", "Male", "Female", "...
## $ HourlyRate               <dbl> 94, 61, 92, 56, 40, 79, 81, 67, 44, 9...
## $ JobInvolvement           <fct> High, Medium, Medium, High, High, Hig...
## $ JobLevel                 <dbl> 2, 2, 1, 1, 1, 1, 1, 1, 3, 2, 1, 2, 1...
## $ JobRole                  <chr> "Sales Executive", "Research Scientis...
## $ JobSatisfaction          <fct> Very High, Medium, High, High, Medium...
## $ MaritalStatus            <chr> "Single", "Married", "Single", "Marri...
## $ MonthlyIncome            <dbl> 5993, 5130, 2090, 2909, 3468, 3068, 2...
## $ MonthlyRate              <dbl> 19479, 24907, 2396, 23159, 16632, 118...
## $ NumCompaniesWorked       <dbl> 8, 1, 6, 1, 9, 0, 4, 1, 0, 6, 0, 0, 1...
## $ Over18                   <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y...
## $ OverTime                 <chr> "Yes", "No", "Yes", "Yes", "No", "No"...
## $ PercentSalaryHike        <dbl> 11, 23, 15, 11, 12, 13, 20, 22, 21, 1...
## $ PerformanceRating        <fct> Excellent, Outstanding, Excellent, Ex...
## $ RelationshipSatisfaction <fct> Low, Very High, Medium, High, Very Hi...
## $ StandardHours            <dbl> 80, 80, 80, 80, 80, 80, 80, 80, 80, 8...
## $ StockOptionLevel         <dbl> 0, 1, 0, 0, 1, 0, 3, 1, 0, 2, 1, 0, 1...
## $ TotalWorkingYears        <dbl> 8, 10, 7, 8, 6, 8, 12, 1, 10, 17, 6, ...
## $ TrainingTimesLastYear    <dbl> 0, 3, 3, 3, 3, 2, 3, 2, 2, 3, 5, 3, 1...
## $ WorkLifeBalance          <fct> Bad, Better, Better, Better, Better, ...
## $ YearsAtCompany           <dbl> 6, 10, 0, 8, 2, 7, 1, 1, 9, 7, 5, 9, ...
## $ YearsInCurrentRole       <dbl> 4, 7, 0, 7, 2, 7, 0, 0, 7, 7, 4, 5, 2...
## $ YearsSinceLastPromotion  <dbl> 0, 1, 0, 3, 2, 3, 0, 0, 1, 7, 0, 0, 4...
## $ YearsWithCurrManager     <dbl> 5, 7, 0, 0, 2, 6, 0, 0, 8, 7, 3, 8, 3...

7.1.3 Factoring Character Data

There are still some character variables:

data_merged_tbl %>% select_if(is.character) %>% glimpse()

## Observations: 1,250
## Variables: 9
## $ Attrition      <chr> "Yes", "No", "Yes", "No", "No", "No", "No", "No...
## $ BusinessTravel <chr> "Travel_Rarely", "Travel_Frequently", "Travel_R...
## $ Department     <chr> "Sales", "Research & Development", "Research & ...
## $ EducationField <chr> "Life Sciences", "Life Sciences", "Other", "Lif...
## $ Gender         <chr> "Female", "Male", "Male", "Female", "Male", "Ma...
## $ JobRole        <chr> "Sales Executive", "Research Scientist", "Labor...
## $ MaritalStatus  <chr> "Single", "Married", "Single", "Married", "Marr...
## $ Over18         <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y...
## $ OverTime       <chr> "Yes", "No", "Yes", "Yes", "No", "No", "Yes", "...

Need to transformed into factors. Evaluate if the features are ordered.

data_merged_tbl %>% distinct(BusinessTravel)

## # A tibble: 3 x 1
##   BusinessTravel   
##   <chr>            
## 1 Travel_Rarely    
## 2 Travel_Frequently
## 3 Non-Travel

The natural order should be: Non_Travel –> Travel_Rarely –> Travel_Frequently. Change it!

Evaluate the levels of the factors:

# Reorder some factors
data_merged_tbl %>% mutate_if(is.character, as.factor) %>%  select_if(is.factor) %>% map(levels) %>% head()

## $Attrition
## [1] "No"  "Yes"
## 
## $BusinessTravel
## [1] "Non-Travel"        "Travel_Frequently" "Travel_Rarely"    
## 
## $Department
## [1] "Human Resources"        "Research & Development"
## [3] "Sales"                 
## 
## $Education
## [1] "Below College" "College"       "Bachelor"      "Master"       
## [5] "Doctor"       
## 
## $EducationField
## [1] "Human Resources"  "Life Sciences"    "Marketing"       
## [4] "Medical"          "Other"            "Technical Degree"
## 
## $EnvironmentSatisfaction
## [1] "Low"       "Medium"    "High"      "Very High"

BusinessTravel and MartitalStatus are the only factor that need to be re-leveled.

data_processed_tbl <- data_merged_tbl %>% mutate_if(is.character, as.factor) %>%
     mutate(BusinessTravel = BusinessTravel %>% fct_relevel("Non-Travel", "Travel_Rarely", "Travel_Frequently"),
            MaritalStatus = MaritalStatus %>% fct_relevel("Single", "Married", "Divorced"))

data_processed_tbl %>% select(BusinessTravel, MaritalStatus) %>% map(levels)

## $BusinessTravel
## [1] "Non-Travel"        "Travel_Rarely"     "Travel_Frequently"
## 
## $MaritalStatus
## [1] "Single"   "Married"  "Divorced"

7.3.1 Custom Function: plot_hist_facet

To inspect feature distributions and identify transformations needed to make them properly formatted for correlation analysis.

• palette()[[3]] is just a hex code “18BC9C”
• as.factor orders factors alphabetically (as_factor maintains the order in which they appear)
• using as.factor and as.numeric covers cases for characters and/or factors - noting special here
• gather does wide df to long df
• fact_rev simply reverses the order of the factors
• facet_wrap splits data into multiple graphs by a categorical columns - or set of columns -. Use functional format to specify facets (~key)

# Plot Faceted Histogram Function ----
plot_hist_facet <- function(data, bins = 10, ncol = 5, fct_reorder = FALSE, fct_rev = FALSE, 
                            fill = palette()[[3]], color = "white", scale = "free"){
     
     data_factored <- data %>% mutate_if(is.character, as.factor) %>% mutate_if(is.factor, as.numeric) %>% 
          gather(key = key, value = value, factor_key = TRUE)
     
     if(fct_reorder){data_factored <- data_factored %>% mutate(key = as.character(key) %>% as.factor())}
     
     if(fct_rev){data_factored <- data_factored %>% mutate(key = fct_rev(key))}
     
     g <- data_factored %>% ggplot(aes(x = value, group = key)) + geom_histogram(bins = bins, 
               fill = fill, color = color) + facet_wrap(~ key, ncol = ncol, scale = scale) + theme_tq()
     return(g)}

Here is the output as an example:

train_raw_tbl %>% select(Attrition, everything()) %>% plot_hist_facet(bins = 10, ncol =5, fct_rev = F)

#ggsave("./images/plot_hist.jpg")

Useful to make the label the first plot

7.3.2 recipes: Preprocessing Data For Machines

The idea of the recipes package is to define a recipe or blueprint that can be used to sequentially define the encodings and preprocessing of the data (i.e. feature engineering).

In recipes, there are no constraints related to the order in which steps are added to the recipe. However, there are some general suggestions that you should consider:

• If using a Box-Cox transformation, don’t center the data first or do any operations that might make the data non-positive. Alternatively, use the Yeo-Johnson transformation so you don’t have to worry about this.
• Recipes do not automatically create dummy variables (unlike most formula methods). If you want to center, scale, or do any other operations on all of the predictors, run step_dummy first so that numeric columns are in the data set instead of factors.
• You should make dummy variables before creating the interactions.
• If you are lumping infrequently categories together with step_other, call step_other before step_dummy.

Here is a suggested order of potential steps that should work for most problems:

1. Impute
2. Individual transformations
   • changes to the data to remove skew (e.g. log), stabilize variance (e.g. Box Cox) or make stationary (e.g. difference for time series)
   • Normality is required for linear models that depend on correlation (correlation analysis, linear and logistic regression)
   • Non-linear models like random forest handle non-linear data because they depend on how the data can be segregated
3. Discretize (if needed and if you have no other choice)
   • discretation can hurt correlations. Often best not to discretize unless there is a specific need to do so
     • Example: Explain the difference between cohorts like Millennials vs Generation X
4. Create dummy variables
5. Create interactions
   • When 2 features have a relationship to each other there are said to interact. Example: The ratio between height and weight of a person is related to obesity. Height and weight interact.
   • Some advanced modeling - deep learning - detect interactions automatically
6. Normalization steps (center, scale, range, etc) Prefer to make this #4 before dummy variables
   • Some ML algorithms (PCA, KNN, Deep Learning) require scaled data
7. Multivariate transformation (e.g. PCA, spatial sign, etc)
   • Useful when data is wide and can be susceptible to over-fitting

For our data, here is the recipe plan we will implement:

• Plan: Correlation Analysis
  • Always do correlation analysis before modeling
  • Impute/Zero Variance
    • Imputation: The act of filling in missing values within features.
      
    • Common methods are:
      • filling in by recency - tidyr: fill
      • filling by similarity - knn impute
  • Transformations
  • Center and scale Dummy Data

We do not have missing values in the current data. However there are features with no variance.

recipes: Setup

First create a recipe from the original data and then specify the processing steps. Recipes can be created manually by sequentially adding roles to variables in a data set. If the analysis only required outcomes and predictors, the easiest way to create the initial recipe is to use the standard formula method:

recipe_obj <- recipe(Attrition ~., data = train_readable_tbl)
recipe_obj

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         34

The data contained in the data argument need not be the training set; this data is only used to catalog the names of the variables and their types (e.g. numeric, etc.).

recipes: Zero Variance Features

To add a step to remove zero variance features, it is simple to do:

recipe_obj %>% step_zv(all_predictors())

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         34
## 
## Operations:
## 
## Zero variance filter on all_predictors()

prep() prepares the recipe performing preliminary calculations. Useful to see what the step is doing.

• prep() does not transform the data - it just determines is is needed to transform
• bake() performs transformation in new data (even if it is the same data you have already used) - bake takes as its newdata argument the name of the dataset typically used in the setup: recipe_obj <- recipe(Attrition ~., data = train_readable_tbl). Useful to apply the recipe to any data set (commonly test and train)
• Essentially 3 steps:
  1. Create the instructions with the recipe and steps
  2. Prepare the recipe (prep())
  3. bake the recipe (`bake(newdata = testdataname)’)

recipe_obj %>% step_zv(all_predictors()) %>% prep()

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         34
## 
## Training data contained 1250 data points and no missing data.
## 
## Operations:
## 
## Zero variance filter removed EmployeeCount, ... [trained]

It is helpful to know 3 features will be removed when the recipe is implemented.

recipes: Transformations

Identify which features have high skewness. Could use the facet plot illustrate above or compute using code (helpful when you have many features):

High Value –> fat tail on left | Low Value –> fat tail on right

skewness is a function from PerformanceAnalytics (package is loaded with tidyquant) computes the skewness of an univariate distribution.

factor_key preserves the order of the variables

train_readable_tbl %>% select_if(is.numeric) %>% map_df(skewness) %>% gather(factor_key = TRUE) %>% arrange(desc(value))

## # A tibble: 19 x 2
##    key                          value
##    <fct>                        <dbl>
##  1 YearsSinceLastPromotion    1.94   
##  2 YearsAtCompany             1.76   
##  3 MonthlyIncome              1.35   
##  4 TotalWorkingYears          1.11   
##  5 NumCompaniesWorked         1.02   
##  6 JobLevel                   1.02   
##  7 DistanceFromHome           0.985  
##  8 StockOptionLevel           0.968  
##  9 YearsInCurrentRole         0.933  
## 10 PercentSalaryHike          0.802  
## 11 YearsWithCurrManager       0.802  
## 12 TrainingTimesLastYear      0.557  
## 13 Age                        0.422  
## 14 EmployeeNumber             0.0121 
## 15 MonthlyRate                0.0100 
## 16 DailyRate                  0.00578
## 17 HourlyRate                -0.0263 
## 18 EmployeeCount            NaN      
## 19 StandardHours            NaN

# NOTE:  high value --> fat tail on left, very low --> fat tail on right

Judgment is used to determine a cutoff for what is/not skewed. Evaluate the values from above:

• Appears to be a big difference between YearsWithCurrManager = 0.802 and TrainingTimesLastYear = 0.55

The facet plot indicates this too. as.character is used so we can plot. . . .

skewed_feature_names <- train_readable_tbl %>% select_if(is.numeric) %>% map_df(skewness) %>% 
     gather(factor_key = TRUE) %>% arrange(desc(value)) %>% filter(value >= 0.8) %>% 
     pull(key) %>% as.character()
skewed_feature_names

##  [1] "YearsSinceLastPromotion" "YearsAtCompany"         
##  [3] "MonthlyIncome"           "TotalWorkingYears"      
##  [5] "NumCompaniesWorked"      "JobLevel"               
##  [7] "DistanceFromHome"        "StockOptionLevel"       
##  [9] "YearsInCurrentRole"      "PercentSalaryHike"      
## [11] "YearsWithCurrManager"

With this information, plot the features:

train_readable_tbl %>% select(skewed_feature_names) %>% plot_hist_facet()

The plot shows JobLevel and StockOptionLevel are factors and therefore are not transformed here (wait for dummying)

factor_names <- c("JobLevel", "StockOptionLevel")

skewed_feature_names <- train_readable_tbl %>% select_if(is.numeric) %>% map_df(skewness) %>% 
     gather(factor_key = TRUE) %>% arrange(desc(value)) %>% filter(value >= 0.8) %>% 
     filter(!key %in% factor_names) %>% pull(key) %>% as.character()

We now have the information we need to continue building our recipe:

recipe_obj <- recipe(Attrition ~., data = train_readable_tbl) %>% step_zv(all_predictors()) %>% 
     step_YeoJohnson(skewed_feature_names) %>% 
     step_num2factor(factor_names)
recipe_obj %>% prep()

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         34
## 
## Training data contained 1250 data points and no missing data.
## 
## Operations:
## 
## Zero variance filter removed EmployeeCount, ... [trained]
## Yeo-Johnson transformation on 9 items [trained]
## Factor variables from JobLevel, StockOptionLevel [trained]

step_YeoJohnson - power transformation. Instead of sqrt, looks for diff roots b/t -5 to 5 to eliminate skew

To make sure the recipe worked, facet plot the numeric features that were the target of the YeoJohnson transformation:

recipe_obj %>% prep() %>% bake(train_readable_tbl) %>% select(skewed_feature_names) %>% plot_hist_facet()

Skewness is significantly reduced!

recipes: center & scale

Required for algos that require scaling including kmeans, Deep Learning, PCA, SVMs.

Scaling does not typically hurt algos - so do it

train_readable_tbl %>% select_if(is.numeric) %>% plot_hist_facet()

Note the x axis values are different from one and other - different ranges. DailyRate would likely dominate!

Must do in order: Center then Scale (C before S)

Note - H2O does not require centering and scaling - more on this later.

recipe_obj <- recipe(Attrition ~., data = train_readable_tbl) %>% step_zv(all_predictors()) %>% 
     step_YeoJohnson(skewed_feature_names) %>% step_num2factor(factor_names) %>% 
     step_center(all_numeric()) %>% 
     step_scale(all_numeric())
recipe_obj

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         34
## 
## Operations:
## 
## Zero variance filter on all_predictors()
## Yeo-Johnson transformation on skewed_feature_names
## Factor variables from factor_names
## Centering for all_numeric()
## Scaling for all_numeric()

The output from the recipe object is helpful but we are interested in understanding what it is really doing behind the scenes:

recipe_obj %>% prep()#Note:  recipe output is a list

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         34
## 
## Training data contained 1250 data points and no missing data.
## 
## Operations:
## 
## Zero variance filter removed EmployeeCount, ... [trained]
## Yeo-Johnson transformation on 9 items [trained]
## Factor variables from JobLevel, StockOptionLevel [trained]
## Centering for Age, DailyRate, ... [trained]
## Scaling for Age, DailyRate, ... [trained]

recipe_obj %>% prep() %>% bake(newdata = train_readable_tbl) %>% select_if(is.numeric) %>% plot_hist_facet()

All the x-axis scales are updated for the independent features!

recipes: dummy variables

One hot encoding - expands categorical variables into separate columns of 1/0. Important to detect patterns in unordered data.

all_nominal() selects only categorical/factor data

To see what output we will get, evaluate only JobRole as our example:

dummied_recipe_obj <- recipe(Attrition ~., data = train_readable_tbl) %>% step_zv(all_predictors()) %>% step_YeoJohnson(skewed_feature_names) %>% 
     step_num2factor(factor_names) %>% step_center(all_numeric()) %>% step_scale(all_numeric()) %>% 
     step_dummy(all_nominal()) %>% prep() %>% bake(newdata = train_readable_tbl) %>% select(contains("JobRole")) %>% plot_hist_facet(ncol = 3)

JobRole has been expanded to 8 columns (1 less than the number of unique values in JobRole. Healthcare Representative Manager was left out).

recipe::step_dummy was updated with a new option - one_hot: A logical. For C levels, should C dummy variables be created rather than C-1? Very useful for plotting!

recipe_obj <- recipe(Attrition ~., data = train_readable_tbl) %>% step_zv(all_predictors()) %>% step_YeoJohnson(skewed_feature_names) %>% 
     step_num2factor(factor_names) %>% step_center(all_numeric()) %>% step_scale(all_numeric()) %>% 
     step_dummy(all_nominal())

recipes: Final Recipe

Transform the data. . .. first review what the recipe will do . . . .

recipe_obj <- recipe_obj %>% prep()
recipe_obj

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         34
## 
## Training data contained 1250 data points and no missing data.
## 
## Operations:
## 
## Zero variance filter removed EmployeeCount, ... [trained]
## Yeo-Johnson transformation on 9 items [trained]
## Factor variables from JobLevel, StockOptionLevel [trained]
## Centering for Age, DailyRate, ... [trained]
## Scaling for Age, DailyRate, ... [trained]
## Dummy variables from BusinessTravel, Department, ... [trained]

train_tbl <- bake(recipe_obj, newdata = train_readable_tbl)

test_tbl <- bake(recipe_obj, newdata = test_readable_tbl)

All of this effort is so we can perform a correlation analysis.

Note that the recipe changes our target label:

class(train_tbl$Attrition_Yes)

## [1] "numeric"

7.4.1 get_cor

A function to take data and measure correlation.

• pairwise.complete.obs is not the default but it is what you nearly always want to use
  • If you use the default - “everything” - risk of getting missing value errors
• fct_reorder reorders factor by another column. Changes the levels of a factor but not the DF - use arrange for that
• correlation analysis only works with numeric data
• quo_name changes the stored value to “Attrition_Yes” from a quosure
• mutate_if(is.character, as.factor) %>% mutate_if(is.factor, as.numeric) prevents function from failing is there is character data
• cor from stats returns a square correlation df correlating every feature against all the others. The return is square meaning the \(number of rows = number of columns\)
  • use is an optional character string giving a method for computing covariances in the presence of missing values.
    • This must be (an abbreviation of) one of the strings “everything”, “all.obs”, “complete.obs”, “na.or.complete”, or “pairwise.complete.obs”
    • If use is “everything”, NAs will propagate conceptually, i.e., a resulting value will be NA whenever one of its contributing observations is NA.
    • If use is “all.obs”, then the presence of missing observations will produce an error. If use is “complete.obs” then missing values are handled by casewise deletion (and if there are no complete cases, that gives an error).
    • “na.or.complete” is the same unless there are no complete cases, that gives NA.
    • If use has the value “pairwise.complete.obs” then the correlation or covariance between each pair of variables is computed using all complete pairs of observations on those variables. This can result in covariance or correlation matrices which are not positive semi-definite, as well as NA entries if there are no complete pairs for that pair of variables.
• The (.) in the names function simply represents the data in the function. This creates a new character column called feature. this makes is square.
• select(feature, !! feature_expr) produces a tibble with 2 columns - the new feature column and the Attrition_Yes values
• filter(!(feature == feature_name)) simply removes the value in the feature column that equals Attrition_Yes - it simply is not needed/redundant
• Recall as_factor maintains the order in which the values appear - not alphabetized
• fct_reorder reorders a factor by another column. In this case, reorder by the Attrition_Yes value. Helpful! Requires arrange() to update the df.
• fct_rev reverses the ordering from the previous line of code

fct_reorder and fct_rev is useful for plotting

get_cor <- function(data, target, use = "pairwise.complete.obs", fct_reorder = FALSE, fct_rev = FALSE){
     
     feature_expr <- enquo(target)
     feature_name <- quo_name(feature_expr)
     
     data_cor <- data %>% mutate_if(is.character, as.factor) %>% mutate_if(is.factor, as.numeric) %>%
          cor(use = use) %>% as.tibble() %>% mutate(feature = names(.)) %>% 
          select(feature, !! feature_expr) %>% filter(!(feature == feature_name)) %>%
          mutate_if(is.character, as_factor)
     
     if(fct_reorder){data_cor <- data_cor %>% mutate(feature = fct_reorder(feature, 
                                                            !! feature_expr)) %>% arrange(feature)}
     
     if(fct_rev){data_cor <- data_cor %>% mutate(feature = fct_rev(feature)) %>% arrange(feature)}
     
     return(data_cor)}

Example of get_cor()

train_tbl %>% get_cor(Attrition_Yes, fct_reorder = T, fct_rev = T)

## Warning in stats::cor(x, ...): the standard deviation is zero

## # A tibble: 65 x 2
##    feature                          Attrition_Yes
##    <fct>                                    <dbl>
##  1 PerformanceRating_Good                 NA     
##  2 OverTime_Yes                            0.240 
##  3 JobRole_Sales.Representative            0.153 
##  4 BusinessTravel_Travel_Frequently        0.103 
##  5 Department_Sales                        0.0891
##  6 DistanceFromHome                        0.0822
##  7 JobRole_Laboratory.Technician           0.0737
##  8 EducationField_Technical.Degree         0.0731
##  9 JobRole_Human.Resources                 0.0718
## 10 JobInvolvement_Medium                   0.0604
## # ... with 55 more rows

• Good correlation with real-world data could be +/-0.1 or greater. It just depends on the problem. It’s very often that the data is messy, but if you get one good feature or a few moderately good features (e.g. 0.2 or more) then this could be good.
• The goal with the preprocessing is to help improve the correlation even if the relationship is nonlinear. For example that is why we do YeoJohnson on some of the skewed data. However, there may be other ways. We could definitely run a random forest and that will pick up any non linearity. We could then try to get feature importance.

7.4.2 Note on Multicolinearity

Multicolinearity is when two features in a model are highly correlated because it relates to whether or not the coefficients are uniquely identified. This relates to explanatory power of the model.

The goals are different for machine learning and traditional statistics, which is why multicolinearity is not traditionally checked for in ML but is in statistics.

Machine Learning: Assessing multicolinearity is not typically a concern in machine learning where the goal is to get a high performance prediction and more features helps in this pursuit. ML models can be linear or nonlinear and typically have regularization which fights over fitting, a situation that can be caused by many features that are correlated. This enables extreme accuracy, but the feature weights (coefficients) get lost, which reduces the explainability of the model.

In machine learning, we will have hundreds of features, but a few of the features will likely dominate for any given observation. This is where LIME comes in. LIME is a tool that allows us to cut through to the core of the features that are driving the explanation for any given model. It approximates nonlinear models locally, which enables us to detect which features are driving the prediction. It’s not impervious to colinear features but this is why we need to understand what features are going into a model from a business-intuition standpoint and discuss features with Subject Matter Experts. See this resource for more discussion on the LIME topic: https://github.com/marcotcr/lime/issues/178

Traditional Statistics: In traditional statistics it’s very common to inspect for multicolinearity and remove features. This is because the goal is to build a model using few features that explain the majority of the variance. When removing features, we can assess multicolinearity

The Debate: Because of the two sides, it’s debatable and it depends on your goals of high performance vs ability to explain with as few features as possible.

My Opinion: When it comes to modeling, my opinion is that we should strive for a high performing model once we know we have good features (features that have some predictive power based on a correlation assessment). When it comes to selecting features, we should meet with Subject Matter Experts and pick reasonable features to build a model around doing our best to not include the same feature twice, which is often where multicolinearity arises from (e.g. including birth year and age in the same model). When it comes to explaining, we should use tools like LIME that have the ability to hone in on the features that are contributing most to the prediction.

For more on this topic: https://stats.stackexchange.com/questions/168622/why-is-multicollinearity-not-checked-in-modern-statistics-machine-learning

7.4.3 plot_cor

The data produced by get_cor() is terrific, but understanding the data quickly is accomplished by creating a visualization.

The function below produces an informative visualization of the get_cor() output. Do not be overwhelmed by the function. Many of the parameters are simply ggplot parameters.

• The beginning of the function looks much like get_cor()
• feature_name_text is simply the Attriiotion_Yes value rounded to the number of digits specified by lbl_precision
• New feature named Correlation is created based on the value of Attrition_Yes and outputs a factor = Positive or Negative
• Experiment with the ggplot code to understand how it works. (First set values for the function parameters.)

plot_cor <- function(data, target, fct_reorder = FALSE, fct_rev = FALSE, include_lbl = TRUE, lbl_precision = 2, lbl_position = "outward",
                     size = 2, line_size = 1, vert_size = 1, color_pos = palette_light()[[1]], color_neg = palette_light()[[2]]) {
     
     feature_expr <- enquo(target)
     feature_name <- quo_name(feature_expr)
     
     data_cor <- data %>%
          get_cor(!! feature_expr, fct_reorder = fct_reorder, fct_rev = fct_rev) %>%
          mutate(feature_name_text = round(!! feature_expr, lbl_precision)) %>%
          mutate(Correlation = case_when((!! feature_expr) >= 0 ~ "Positive", TRUE ~ "Negative") %>% as.factor())
     
     g <- data_cor %>% ggplot(aes_string(x = feature_name, y = "feature", group = "feature")) +
          geom_point(aes(color = Correlation), size = size) +
          geom_segment(aes(xend = 0, yend = feature, color = Correlation), size = line_size) +
          geom_vline(xintercept = 0, color = palette_light()[[1]], size = vert_size) +
          expand_limits(x = c(-1, 1)) + theme_tq() + scale_color_manual(values = c(color_neg, color_pos)) 
     
     if (include_lbl) g <- g + geom_label(aes(label = feature_name_text), hjust = lbl_position)
     
     return(g)}

train_tbl %>% plot_cor(target = Attrition_Yes, fct_reorder = T, fct_rev = F)

#ggsave("./images/plot_cor.jpg")

7.4.4 Evaluating plot_cor Output

Data have been prepared and we developed a visualization function to easily evaluate feature correlations.

1. Descriptive Features: age, gender, marital status, etc

train_tbl %>% select(Attrition_Yes, Age, contains("Gender"), contains("MaritalStatus"), NumCompaniesWorked, contains("Over18"), DistanceFromHome) %>% 
     plot_cor(target = Attrition_Yes, fct_reorder = T, fct_rev = F)

Observations:

• As age increases, less likely to leave company
  • DistanceFromHome has a relationship as dsitance from home supports someone more likely to leaeave the company

2. Employment Features

train_tbl %>% select(Attrition_Yes, contains("employee"), contains("department"), contains("job")) %>% 
plot_cor(target = Attrition_Yes, fct_reorder = T, fct_rev = F)

# contains is NOT case sentitive

Observations:

• Certain jobs either support or contradict the lieklihood to leave

3. Compensation Features

train_tbl %>% select(Attrition_Yes, contains("income"), contains("rate"), contains("salary"), 
                     contains("stock")) %>% 
     plot_cor(target = Attrition_Yes, fct_reorder = T, fct_rev = F)

Observations:

-MonhtlyIncome StockOptionLevel are impactful

4. Survey Results

train_tbl %>% select(Attrition_Yes, contains("satisfaction"), contains("life")) %>% 
     plot_cor(target = Attrition_Yes, fct_reorder = T, fct_rev = F)

5. Performance Data

train_tbl %>% select(Attrition_Yes, contains("performance"), contains("involvement")) %>% 
     plot_cor(target = Attrition_Yes, fct_reorder = T, fct_rev = F)

6. Work-Life Features

train_tbl %>% select(Attrition_Yes, contains("overtime"), contains("travel")) %>% 
     plot_cor(target = Attrition_Yes, fct_reorder = T, fct_rev = F)

7. Training and Education

train_tbl %>% select(Attrition_Yes, contains("training"), contains("education")) %>% 
     plot_cor(target = Attrition_Yes, fct_reorder = T, fct_rev = F)

8. Time-Based Features

train_tbl %>% select(Attrition_Yes, contains("years")) %>% plot_cor(target = Attrition_Yes, 
                                                            fct_reorder = T, fct_rev = F)

8.1.1 Inspecting The Leaderboard

There is much data produced by the individual models - the leader is shown below:

automl_models_h2o@leader

## Model Details:
## ==============
## 
## H2OBinomialModel: glm
## Model ID:  GLM_grid_0_AutoML_20181028_212737_model_0 
## GLM Model: summary
##     family  link            regularization
## 1 binomial logit Ridge ( lambda = 0.0042 )
##                                                                  lambda_search
## 1 nlambda = 30, lambda.max = 6.2463, lambda.min = 0.0042, lambda.1se = 0.02055
##   number_of_predictors_total number_of_active_predictors
## 1                         79                          79
##   number_of_iterations                    training_frame
## 1                   54 _804567d077a2a0701a19fb84ffc740ad
## 
## Coefficients: glm coefficients
##                               names coefficients standardized_coefficients
## 1                         Intercept     0.849792                 -1.417281
## 2 JobRole.Healthcare Representative     0.066056                  0.066056
## 3           JobRole.Human Resources     0.167535                  0.167535
## 4     JobRole.Laboratory Technician     0.294450                  0.294450
## 5                   JobRole.Manager    -0.188477                 -0.188477
## 
## ---
##                      names coefficients standardized_coefficients
## 75       TotalWorkingYears    -0.042385                 -0.330973
## 76   TrainingTimesLastYear    -0.180422                 -0.231657
## 77          YearsAtCompany     0.102775                  0.649596
## 78      YearsInCurrentRole    -0.107887                 -0.399926
## 79 YearsSinceLastPromotion     0.098246                  0.321247
## 80    YearsWithCurrManager    -0.118166                 -0.425507
## 
## H2OBinomialMetrics: glm
## ** Reported on training data. **
## 
## MSE:  0.077
## RMSE:  0.28
## LogLoss:  0.27
## Mean Per-Class Error:  0.19
## AUC:  0.9
## Gini:  0.8
## R^2:  0.44
## Residual Deviance:  573
## AIC:  733
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error       Rate
## No     838  51 0.057368    =51/889
## Yes     58 118 0.329545    =58/176
## Totals 896 169 0.102347  =109/1065
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.348528 0.684058 128
## 2                       max f2  0.173866 0.739300 220
## 3                 max f0point5  0.517226 0.767857  79
## 4                 max accuracy  0.517226 0.906103  79
## 5                max precision  0.962484 1.000000   0
## 6                   max recall  0.006173 1.000000 388
## 7              max specificity  0.962484 1.000000   0
## 8             max absolute_mcc  0.348528 0.623206 128
## 9   max min_per_class_accuracy  0.189074 0.825647 210
## 10 max mean_per_class_accuracy  0.218153 0.836509 190
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## H2OBinomialMetrics: glm
## ** Reported on validation data. **
## 
## MSE:  0.086
## RMSE:  0.29
## LogLoss:  0.3
## Mean Per-Class Error:  0.24
## AUC:  0.81
## Gini:  0.62
## R^2:  0.26
## Residual Deviance:  112
## AIC:  272
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error     Rate
## No     148  12 0.075000  =12/160
## Yes     10  15 0.400000   =10/25
## Totals 158  27 0.118919  =22/185
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.314235 0.576923  26
## 2                       max f2  0.167744 0.620915  52
## 3                 max f0point5  0.648042 0.612245   5
## 4                 max accuracy  0.648042 0.897297   5
## 5                max precision  0.933029 1.000000   0
## 6                   max recall  0.010212 1.000000 158
## 7              max specificity  0.933029 1.000000   0
## 8             max absolute_mcc  0.314235 0.508369  26
## 9   max min_per_class_accuracy  0.167744 0.760000  52
## 10 max mean_per_class_accuracy  0.167744 0.773750  52
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## H2OBinomialMetrics: glm
## ** Reported on cross-validation data. **
## ** 5-fold cross-validation on training data (Metrics computed for combined holdout predictions) **
## 
## MSE:  0.094
## RMSE:  0.31
## LogLoss:  0.32
## Mean Per-Class Error:  0.24
## AUC:  0.85
## Gini:  0.69
## R^2:  0.32
## Residual Deviance:  684
## AIC:  844
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error       Rate
## No     799  90 0.101237    =90/889
## Yes     65 111 0.369318    =65/176
## Totals 864 201 0.145540  =155/1065
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.303861 0.588859 155
## 2                       max f2  0.192745 0.659068 208
## 3                 max f0point5  0.481045 0.651408  86
## 4                 max accuracy  0.481045 0.881690  86
## 5                max precision  0.949434 1.000000   0
## 6                   max recall  0.002279 1.000000 396
## 7              max specificity  0.949434 1.000000   0
## 8             max absolute_mcc  0.359402 0.511298 127
## 9   max min_per_class_accuracy  0.173231 0.772727 219
## 10 max mean_per_class_accuracy  0.192745 0.781103 208
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## Cross-Validation Metrics Summary: 
##                 mean          sd cv_1_valid cv_2_valid cv_3_valid
## accuracy  0.87323946  0.02347418 0.84976524 0.82159626 0.90140843
## auc        0.8425771 0.026760576  0.8296127  0.8300225  0.9127907
## err       0.12676056  0.02347418 0.15023474 0.17840375 0.09859155
## err_count       27.0         5.0       32.0       38.0       21.0
## f0point5  0.63031274  0.07437051  0.6377551 0.45064378  0.7763975
##            cv_4_valid cv_5_valid
## accuracy   0.91079813  0.8826291
## auc         0.7989362 0.84152335
## err       0.089201875 0.11737089
## err_count        19.0       25.0
## f0point5    0.6179775  0.6687898
## 
## ---
##                         mean          sd cv_1_valid cv_2_valid cv_3_valid
## precision          0.6580986 0.096957624 0.65789473  0.4117647  0.8333333
## r2                0.30837095 0.061606336  0.2931284 0.25788376  0.4637525
## recall             0.5819287  0.06444905  0.5681818  0.7241379  0.6097561
## residual_deviance  135.58003   11.968308  165.20996  127.08945  117.54124
## rmse              0.30468094 0.014051158 0.34037712  0.2954366 0.28870824
## specificity        0.9306463 0.035521794  0.9230769  0.8369565  0.9709302
##                   cv_4_valid cv_5_valid
## precision             0.6875        0.7
## r2                0.20434217 0.32274798
## recall                  0.44  0.5675676
## residual_deviance   125.3069  142.75262
## rmse              0.28709954 0.31178322
## specificity        0.9734042  0.9488636

H₂O even provides a leaderboard - how cool is that?

automl_models_h2o@leaderboard

##                                                model_id  auc logloss
## 1             GLM_grid_0_AutoML_20181028_212737_model_0 0.85    0.32
## 2 StackedEnsemble_BestOfFamily_0_AutoML_20181028_212737 0.84    0.32
## 3    StackedEnsemble_AllModels_0_AutoML_20181028_212737 0.84    0.32
## 4             GBM_grid_0_AutoML_20181028_212737_model_4 0.82    0.34
## 5             GBM_grid_0_AutoML_20181028_212737_model_5 0.80    0.36
## 6             GBM_grid_0_AutoML_20181028_212737_model_1 0.80    0.36
##   mean_per_class_error rmse   mse
## 1                 0.24 0.31 0.094
## 2                 0.23 0.31 0.093
## 3                 0.24 0.30 0.093
## 4                 0.25 0.32 0.103
## 5                 0.28 0.33 0.109
## 6                 0.27 0.33 0.107
## 
## [21 rows x 6 columns]

8.1.2 Extracting Models From The Leaderboard

Getting the leader is pretty simple:

automl_models_h2o@leader

## Model Details:
## ==============
## 
## H2OBinomialModel: glm
## Model ID:  GLM_grid_0_AutoML_20181028_212737_model_0 
## GLM Model: summary
##     family  link            regularization
## 1 binomial logit Ridge ( lambda = 0.0042 )
##                                                                  lambda_search
## 1 nlambda = 30, lambda.max = 6.2463, lambda.min = 0.0042, lambda.1se = 0.02055
##   number_of_predictors_total number_of_active_predictors
## 1                         79                          79
##   number_of_iterations                    training_frame
## 1                   54 _804567d077a2a0701a19fb84ffc740ad
## 
## Coefficients: glm coefficients
##                               names coefficients standardized_coefficients
## 1                         Intercept     0.849792                 -1.417281
## 2 JobRole.Healthcare Representative     0.066056                  0.066056
## 3           JobRole.Human Resources     0.167535                  0.167535
## 4     JobRole.Laboratory Technician     0.294450                  0.294450
## 5                   JobRole.Manager    -0.188477                 -0.188477
## 
## ---
##                      names coefficients standardized_coefficients
## 75       TotalWorkingYears    -0.042385                 -0.330973
## 76   TrainingTimesLastYear    -0.180422                 -0.231657
## 77          YearsAtCompany     0.102775                  0.649596
## 78      YearsInCurrentRole    -0.107887                 -0.399926
## 79 YearsSinceLastPromotion     0.098246                  0.321247
## 80    YearsWithCurrManager    -0.118166                 -0.425507
## 
## H2OBinomialMetrics: glm
## ** Reported on training data. **
## 
## MSE:  0.077
## RMSE:  0.28
## LogLoss:  0.27
## Mean Per-Class Error:  0.19
## AUC:  0.9
## Gini:  0.8
## R^2:  0.44
## Residual Deviance:  573
## AIC:  733
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error       Rate
## No     838  51 0.057368    =51/889
## Yes     58 118 0.329545    =58/176
## Totals 896 169 0.102347  =109/1065
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.348528 0.684058 128
## 2                       max f2  0.173866 0.739300 220
## 3                 max f0point5  0.517226 0.767857  79
## 4                 max accuracy  0.517226 0.906103  79
## 5                max precision  0.962484 1.000000   0
## 6                   max recall  0.006173 1.000000 388
## 7              max specificity  0.962484 1.000000   0
## 8             max absolute_mcc  0.348528 0.623206 128
## 9   max min_per_class_accuracy  0.189074 0.825647 210
## 10 max mean_per_class_accuracy  0.218153 0.836509 190
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## H2OBinomialMetrics: glm
## ** Reported on validation data. **
## 
## MSE:  0.086
## RMSE:  0.29
## LogLoss:  0.3
## Mean Per-Class Error:  0.24
## AUC:  0.81
## Gini:  0.62
## R^2:  0.26
## Residual Deviance:  112
## AIC:  272
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error     Rate
## No     148  12 0.075000  =12/160
## Yes     10  15 0.400000   =10/25
## Totals 158  27 0.118919  =22/185
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.314235 0.576923  26
## 2                       max f2  0.167744 0.620915  52
## 3                 max f0point5  0.648042 0.612245   5
## 4                 max accuracy  0.648042 0.897297   5
## 5                max precision  0.933029 1.000000   0
## 6                   max recall  0.010212 1.000000 158
## 7              max specificity  0.933029 1.000000   0
## 8             max absolute_mcc  0.314235 0.508369  26
## 9   max min_per_class_accuracy  0.167744 0.760000  52
## 10 max mean_per_class_accuracy  0.167744 0.773750  52
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## H2OBinomialMetrics: glm
## ** Reported on cross-validation data. **
## ** 5-fold cross-validation on training data (Metrics computed for combined holdout predictions) **
## 
## MSE:  0.094
## RMSE:  0.31
## LogLoss:  0.32
## Mean Per-Class Error:  0.24
## AUC:  0.85
## Gini:  0.69
## R^2:  0.32
## Residual Deviance:  684
## AIC:  844
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error       Rate
## No     799  90 0.101237    =90/889
## Yes     65 111 0.369318    =65/176
## Totals 864 201 0.145540  =155/1065
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.303861 0.588859 155
## 2                       max f2  0.192745 0.659068 208
## 3                 max f0point5  0.481045 0.651408  86
## 4                 max accuracy  0.481045 0.881690  86
## 5                max precision  0.949434 1.000000   0
## 6                   max recall  0.002279 1.000000 396
## 7              max specificity  0.949434 1.000000   0
## 8             max absolute_mcc  0.359402 0.511298 127
## 9   max min_per_class_accuracy  0.173231 0.772727 219
## 10 max mean_per_class_accuracy  0.192745 0.781103 208
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## Cross-Validation Metrics Summary: 
##                 mean          sd cv_1_valid cv_2_valid cv_3_valid
## accuracy  0.87323946  0.02347418 0.84976524 0.82159626 0.90140843
## auc        0.8425771 0.026760576  0.8296127  0.8300225  0.9127907
## err       0.12676056  0.02347418 0.15023474 0.17840375 0.09859155
## err_count       27.0         5.0       32.0       38.0       21.0
## f0point5  0.63031274  0.07437051  0.6377551 0.45064378  0.7763975
##            cv_4_valid cv_5_valid
## accuracy   0.91079813  0.8826291
## auc         0.7989362 0.84152335
## err       0.089201875 0.11737089
## err_count        19.0       25.0
## f0point5    0.6179775  0.6687898
## 
## ---
##                         mean          sd cv_1_valid cv_2_valid cv_3_valid
## precision          0.6580986 0.096957624 0.65789473  0.4117647  0.8333333
## r2                0.30837095 0.061606336  0.2931284 0.25788376  0.4637525
## recall             0.5819287  0.06444905  0.5681818  0.7241379  0.6097561
## residual_deviance  135.58003   11.968308  165.20996  127.08945  117.54124
## rmse              0.30468094 0.014051158 0.34037712  0.2954366 0.28870824
## specificity        0.9306463 0.035521794  0.9230769  0.8369565  0.9709302
##                   cv_4_valid cv_5_valid
## precision             0.6875        0.7
## r2                0.20434217 0.32274798
## recall                  0.44  0.5675676
## residual_deviance   125.3069  142.75262
## rmse              0.28709954 0.31178322
## specificity        0.9734042  0.9488636

To evaluate a model other than the leader, you must use the model_id in the leaderboard:

#h2o.getModel("StackedEnsemble_BestOfFamily_0_AutoML_20181028_112225")

The code above is commented only because when the RMD is knitted, different models will be built and the model_id will change.

8.1.3 Custom Function: extract_h2o_model_by_position

extract_h2o_model_name_by_position <- function(h2o_leaderboard, position = 1, verbose = TRUE){
     
     model_name <- h2o_leaderboard %>% as.tibble() %>% slice(position) %>% pull(model_id)
     if(verbose) message(model_name)
     return(model_name)}

# Test new function
extract_h2o_model_name_by_position(automl_models_h2o@leaderboard, position = 3, verbose = T) %>% h2o.getModel()

## StackedEnsemble_AllModels_0_AutoML_20181028_212737

## Model Details:
## ==============
## 
## H2OBinomialModel: stackedensemble
## Model ID:  StackedEnsemble_AllModels_0_AutoML_20181028_212737 
## NULL
## 
## 
## H2OBinomialMetrics: stackedensemble
## ** Reported on training data. **
## 
## MSE:  0.051
## RMSE:  0.23
## LogLoss:  0.19
## Mean Per-Class Error:  0.13
## AUC:  0.96
## Gini:  0.93
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error      Rate
## No     863  26 0.029246   =26/889
## Yes     39 137 0.221591   =39/176
## Totals 902 163 0.061033  =65/1065
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.372597 0.808260 140
## 2                       max f2  0.230229 0.826039 179
## 3                 max f0point5  0.621124 0.860927  90
## 4                 max accuracy  0.419274 0.938967 132
## 5                max precision  0.983435 1.000000   0
## 6                   max recall  0.058457 1.000000 336
## 7              max specificity  0.983435 1.000000   0
## 8             max absolute_mcc  0.372597 0.772834 140
## 9   max min_per_class_accuracy  0.168999 0.889764 210
## 10 max mean_per_class_accuracy  0.230229 0.895794 179
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## H2OBinomialMetrics: stackedensemble
## ** Reported on validation data. **
## 
## MSE:  0.083
## RMSE:  0.29
## LogLoss:  0.29
## Mean Per-Class Error:  0.21
## AUC:  0.82
## Gini:  0.65
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error     Rate
## No     143  17 0.106250  =17/160
## Yes      8  17 0.320000    =8/25
## Totals 151  34 0.135135  =25/185
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.225906 0.576271  33
## 2                       max f2  0.092863 0.644172  62
## 3                 max f0point5  0.850666 0.612245   5
## 4                 max accuracy  0.850666 0.897297   5
## 5                max precision  0.982857 1.000000   0
## 6                   max recall  0.041520 1.000000 146
## 7              max specificity  0.982857 1.000000   0
## 8             max absolute_mcc  0.225906 0.506436  33
## 9   max min_per_class_accuracy  0.105688 0.760000  55
## 10 max mean_per_class_accuracy  0.092863 0.788750  62
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## H2OBinomialMetrics: stackedensemble
## ** Reported on cross-validation data. **
## ** 5-fold cross-validation on training data (Metrics computed for combined holdout predictions) **
## 
## MSE:  0.093
## RMSE:  0.3
## LogLoss:  0.32
## Mean Per-Class Error:  0.24
## AUC:  0.84
## Gini:  0.68
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##         No Yes    Error       Rate
## No     807  82 0.092238    =82/889
## Yes     67 109 0.380682    =67/176
## Totals 874 191 0.139906  =149/1065
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.278936 0.594005 157
## 2                       max f2  0.150185 0.645325 219
## 3                 max f0point5  0.505251 0.667808  88
## 4                 max accuracy  0.505251 0.885446  88
## 5                max precision  0.982169 1.000000   0
## 6                   max recall  0.033742 1.000000 393
## 7              max specificity  0.982169 1.000000   0
## 8             max absolute_mcc  0.437944 0.526086 102
## 9   max min_per_class_accuracy  0.113342 0.755682 253
## 10 max mean_per_class_accuracy  0.150185 0.774744 219
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`

8.1.4 Saving & Loading H2O Models

Saving models is simple. First use h2o.getModel then save it by running h2o.saveModel. Easy! Below the code is not evaluated because each time the RMD is knitted, the models will change.

h2o.getModel("GLM_grid_0_AutoML_20181028_112225_model_0") %>%
     h2o.saveModel(path = "04_Modeling/h2o_models/")

h2o.getModel("StackedEnsemble_BestOfFamily_0_AutoML_20181028_112225") %>%
     h2o.saveModel(path = "04_Modeling/h2o_models/")

Loaded as saved model is equally simple using h2o.loadModel:

#h2o.loadModel("04_Modeling/h2o_models/GLM_grid_0_AutoML_20181028_112225_model_0")

getModel returns a reference to an existing model in the H2O instance. loadModel loads a saved H2O model from disk.

8.1.5 Making Predictions

Load a model:

stacked_ensemble_h2o <- h2o.loadModel("04_Modeling/h2o_models/GLM_grid_0_AutoML_20181028_112225_model_0")

Make predictions with the model using test_tbl and h2o.predict - again, simple!

predictions <- h2o.predict(stacked_ensemble_h2o, newdata = as.h2o(test_tbl))
predictions_tbl <- predictions %>% as.tibble() 
predictions_tbl %>% head()

## # A tibble: 6 x 3
##   predict    No      Yes
##   <fct>   <dbl>    <dbl>
## 1 No      0.981 0.0195  
## 2 No      0.803 0.197   
## 3 No      0.999 0.000884
## 4 No      0.974 0.0260  
## 5 Yes     0.338 0.662   
## 6 No      0.923 0.0766

For H₂O Binary Classification, H₂O provides 3 columns of data:

1. Class prediction
2. 1st class probability
3. 2nd class probability

8.2.1 Train, Validation, & Leaderboard Frames

The models are trained on the training data and then sends the results to the validation data set and learns from the differences.

Repeated Note: The leaderboard_frame is a legacy argument that is not necessary and the data is better served being used for the training and/or validation. Therefore, you should only perform one split (training and validation) as opposed to two (training, validation, test - used for leaderboard_frame).

8.2.2 H2O AutoML Model Parameters

When extracting the model parameters, there sufficient information so the model can be run outside of the H2o framework.H2O automated ML model leaderboard returns a bunch of models. Each model can be extracted and the contents can be opened up to determine the exact model parameters. The only one that will be difficult to recreate is stacked ensemble since this is a combination of models. Others such as GLM, GBM, and Deep Learning can be inspected and recreated. The Bonus Lecture on Grid Search.

automl_models_h2o@leader@parameters

## $model_id
## [1] "GLM_grid_0_AutoML_20181028_212737_model_0"
## 
## $training_frame
## [1] "_804567d077a2a0701a19fb84ffc740ad"
## 
## $validation_frame
## [1] "RTMP_sid_bdf4_5"
## 
## $nfolds
## [1] 5
## 
## $seed
## [1] -5728624559962048512
## 
## $keep_cross_validation_predictions
## [1] TRUE
## 
## $fold_assignment
## [1] "Modulo"
## 
## $family
## [1] "binomial"
## 
## $solver
## [1] "COORDINATE_DESCENT"
## 
## $alpha
## [1] 0.0 0.2 0.4 0.6 0.8 1.0
## 
## $lambda
##  [1] 6.2463 4.5467 3.3095 2.4090 1.7535 1.2763 0.9290 0.6762 0.4922 0.3583
## [11] 0.2608 0.1898 0.1382 0.1006 0.0732 0.0533 0.0388 0.0282 0.0206 0.0150
## [21] 0.0109 0.0079 0.0058 0.0042 0.0031 0.0022 0.0016
## 
## $lambda_search
## [1] TRUE
## 
## $nlambdas
## [1] 30
## 
## $max_iterations
## [1] 300
## 
## $objective_epsilon
## [1] 0.0001
## 
## $gradient_epsilon
## [1] 0.000001
## 
## $link
## [1] "logit"
## 
## $lambda_min_ratio
## [1] 0.0001
## 
## $max_active_predictors
## [1] 5000
## 
## $obj_reg
## [1] 0.00094
## 
## $max_runtime_secs
## [1] 4.3
## 
## $x
##  [1] "JobRole"                  "EducationField"          
##  [3] "JobLevel"                 "Education"               
##  [5] "EnvironmentSatisfaction"  "JobInvolvement"          
##  [7] "JobSatisfaction"          "RelationshipSatisfaction"
##  [9] "StockOptionLevel"         "WorkLifeBalance"         
## [11] "Department"               "BusinessTravel"          
## [13] "MaritalStatus"            "Gender"                  
## [15] "OverTime"                 "PerformanceRating"       
## [17] "Age"                      "DailyRate"               
## [19] "DistanceFromHome"         "EmployeeNumber"          
## [21] "HourlyRate"               "MonthlyIncome"           
## [23] "MonthlyRate"              "NumCompaniesWorked"      
## [25] "PercentSalaryHike"        "TotalWorkingYears"       
## [27] "TrainingTimesLastYear"    "YearsAtCompany"          
## [29] "YearsInCurrentRole"       "YearsSinceLastPromotion" 
## [31] "YearsWithCurrManager"    
## 
## $y
## [1] "Attrition"

automl_models_h2o@leader@parameters$nfolds

## [1] 5

8.2.3 Cross Validation (K-Fold CV)

H₂O Documentation

nfolds is K-fold validation. K-fold validation is used in model development to determine what parameters affect prediction performance by performing multiple sub-tests rather than just one test. The end result is parameter selections based on model performance stability that generalizes better with new data.

For example, for nfolds=5, 6 models are built. The first 5 models (cross-validation models) are built on 80% of the training data, and a different 20% is held out for each of the 5 models. Then the main model is built on 100% of the training data. This main model is the model you get back from H2O in R, Python and Flow (though the CV models are also stored and available to access later).

All 5 cross-validation models contain training metrics (from the 80% training data) and validation metrics (from their 20% holdout/validation data). To compute their individual validation metrics, each of the 5 cross-validation models had to make predictions on their 20% of of rows of the original training frame, and score against the true labels of the 20% holdout.

For the main model, this is how the cross-validation metrics are computed: The 5 holdout predictions are combined into one prediction for the full training dataset (i.e., predictions for every row of the training data, but the model making the prediction for a particular row has not seen that row during training). This “holdout prediction” is then scored against the true labels, and the overall cross-validation metrics are computed.

Each cv-model produces a prediction frame pertaining to its fold. It can be saved and probed from the various clients if keep_cross_validation_predictions parameter is set in the model constructor.

These holdout predictions have some interesting properties. First they have names like:

prediction_GBM_model_1452035702801_1_cv_1

and they contain, unsurprisingly, predictions for the data held out in the fold. They also have the same number of rows as the entire input training frame with 0s filled in for all rows that are not in the hold out.

Retrieve the AUC for a classifier with the xval argument to retrieve the average cross validation AUC:

h2o.auc(automl_models_h2o@leader, train = T, valid = T, xval = T)# xval-cross-validation AUC

## train valid  xval 
##  0.90  0.81  0.85

8.2.4 Grid/Hyperparameter Search

H₂O Documetnation

H2O supports two types of grid search – traditional (or “cartesian”) grid search and random grid search. In a cartesian grid search, users specify a set of values for each hyperparameter that they want to search over, and H2O will train a model for every combination of the hyperparameter values. This means that if you have three hyperparameters and you specify 5, 10 and 2 values for each, your grid will contain a total of 5102 = 100 models.

In random grid search, the user specifies the hyperparameter space in the exact same way, except H2O will sample uniformly from the set of all possible hyperparameter value combinations. In random grid search, the user also specifies a stopping criterion, which controls when the random grid search is completed. The user can tell the random grid search to stop by specifying a maximum number of models or the maximum number of seconds allowed for the search. The user may also specify a performance-metric-based stopping criterion, which will stop the random grid search when the performance stops improving by a specified amount.

Once the grid search is complete, the user can query the grid object and sort the models by a particular performance metric (for example, “AUC”). All models are stored in the H2O cluster and are accessible by model id.

Grid Search

Note the use of h2o.grid and hyper_params

# GBM hyperparamters
gbm_params1 <- list(learn_rate = c(0.01, 0.1), max_depth = c(3, 5, 9), sample_rate = c(0.8, 1.0),
                    col_sample_rate = c(0.2, 0.5, 1.0))

# Train and validate a cartesian grid of GBMs
gbm_grid1 <- h2o.grid("gbm", x = x, y = y, grid_id = "gbm_grid1", training_frame = train, 
                      validation_frame = valid, ntrees = 100, seed = 1,
                      hyper_params = gbm_params1)

# Get the grid results, sorted by validation AUC
gbm_gridperf1 <- h2o.getGrid(grid_id = "gbm_grid1", sort_by = "auc", decreasing = TRUE)
print(gbm_gridperf1)

# Grab the top GBM model, chosen by validation AUC
best_gbm1 <- h2o.getModel(gbm_gridperf1@model_ids[[1]])

# Now let's evaluate the model performance on a test set
# so we get an honest estimate of top model performance
best_gbm_perf1 <- h2o.performance(model = best_gbm1,
                                  newdata = test)
h2o.auc(best_gbm_perf1)  # 0.7781932

# Look at the hyperparamters for the best model
print(best_gbm1@model[["model_summary"]])

Random Grid Search

Note the use of seq, search_criteria: search_criteria <- list(strategy = "RandomDiscrete", max_models = 36, seed = 1) and hyper_parameters

gbm_params2 <- list(learn_rate = seq(0.01, 0.1, 0.01), ax_depth = seq(2, 10, 1),
                    sample_rate = seq(0.5, 1.0, 0.1), col_sample_rate = seq(0.1, 1.0, 0.1))
search_criteria <- list(strategy = "RandomDiscrete", max_models = 36, seed = 1)

# Train and validate a random grid of GBMs
gbm_grid2 <- h2o.grid("gbm", x = x, y = y, grid_id = "gbm_grid2", training_frame = train,
                      validation_frame = valid, ntrees = 100, eed = 1,
                      hyper_params = gbm_params2, earch_criteria = search_criteria)

gbm_gridperf2 <- h2o.getGrid(grid_id = "gbm_grid2", srt_by = "auc", decreasing = TRUE)
print(gbm_gridperf2)

# Grab the top GBM model, chosen by validation AUC
best_gbm2 <- h2o.getModel(gbm_gridperf2@model_ids[[1]])

# Now let's evaluate the model performance on a test set
# so we get an honest estimate of top model performance
best_gbm_perf2 <- h2o.performance(model = best_gbm2,
                                  newdata = test)
h2o.auc(best_gbm_perf2)  # 0.7811332

# Look at the hyperparamters for the best model
print(best_gbm2@model[["model_summary"]])

Specify the RandomDiscrete strategy to perform a random search of all the combinations of your hyperparameters. RandomDiscrete should be combined with at least one early stopping criterion, max_models and/or max_runtime_secs

Examples:

list(strategy = "RandomDiscrete", max_models = 10, seed = 1)
list(strategy = "RandomDiscrete", max_runtime_secs = 3600)
list(strategy = "RandomDiscrete", max_models = 42, max_runtime_secs = 28800)
list(strategy = "RandomDiscrete", stopping_tolerance = 0.001, stopping_rounds = 10)
list(strategy = "RandomDiscrete", 
     stopping_metric = "misclassification", stopping_tolerance = 0.0005, stopping_rounds = 5)

To learn everything about this topic, must visit this tutorial

H₂O deep learning model sets nfolds = 0 by default!

8.3.1 Custom Function: plot_h2o_leaderboard

• str_split simplify makes the output a matrix rather that a list. Makes it easier to extract data [,1]. This allows us to get the type of model from each model ID in the leaderboard.
• rownames_to_column adds row names of a data frame to a column. Gives us a column with the row name as a character. In this case it is very help because it identifies the order in which the models appear in the leaderboard.
• Note the use of as_factor and as.factor. 'as_factor allows us to order the data by AUC. Recall as.factor would order alphabetically.
• Use gather to take advantage of ggplot color and group aesthetics. Data must be in long format. Long format has each faceting values stacked on top of each other making the data frame long rather than wide. - the key column will have the valuesaucandlogloss. Thevaluecolumn will have the values ofaucandlogloss. gather transposes auc and logloss from two columns to one - key and their associated values will be in the value column - Before:

automl_models_h2o@leaderboard %>% as.tibble() %>% 
     mutate(model_type = str_split(model_id, "_", simplify = T)[,1]) %>% slice(1:10) %>% 
     rownames_to_column() %>% 
     mutate(model_id = as_factor(model_id) %>% reorder(auc), model_type = as.factor(model_type)) 

## # A tibble: 10 x 8
##    rowname model_id   auc logloss mean_per_class_~  rmse    mse model_type
##    <chr>   <fct>    <dbl>   <dbl>            <dbl> <dbl>  <dbl> <fct>     
##  1 1       GLM_gri~ 0.846   0.321            0.235 0.307 0.0940 GLM       
##  2 2       Stacked~ 0.840   0.316            0.232 0.306 0.0934 StackedEn~
##  3 3       Stacked~ 0.839   0.316            0.236 0.305 0.0929 StackedEn~
##  4 4       GBM_gri~ 0.816   0.345            0.254 0.320 0.103  GBM       
##  5 5       GBM_gri~ 0.800   0.359            0.285 0.330 0.109  GBM       
##  6 6       GBM_gri~ 0.797   0.362            0.273 0.327 0.107  GBM       
##  7 7       GBM_gri~ 0.795   0.351            0.259 0.322 0.104  GBM       
##  8 8       GBM_gri~ 0.783   0.372            0.293 0.332 0.110  GBM       
##  9 9       GBM_gri~ 0.783   0.379            0.288 0.337 0.113  GBM       
## 10 10      GBM_gri~ 0.781   0.405            0.297 0.353 0.125  GBM

• After:

automl_models_h2o@leaderboard %>% as.tibble() %>% 
     mutate(model_type = str_split(model_id, "_", simplify = T)[,1]) %>% slice(1:10) %>% 
     rownames_to_column() %>% 
     mutate(model_id = as_factor(model_id) %>% reorder(auc), model_type = as.factor(model_type)) %>% 
     gather(key = key, value = value, -c(model_id, model_type, rowname), factor_key = T)

## # A tibble: 50 x 5
##    rowname model_id                              model_type    key   value
##    <chr>   <fct>                                 <fct>         <fct> <dbl>
##  1 1       GLM_grid_0_AutoML_20181028_212737_mo~ GLM           auc   0.846
##  2 2       StackedEnsemble_BestOfFamily_0_AutoM~ StackedEnsem~ auc   0.840
##  3 3       StackedEnsemble_AllModels_0_AutoML_2~ StackedEnsem~ auc   0.839
##  4 4       GBM_grid_0_AutoML_20181028_212737_mo~ GBM           auc   0.816
##  5 5       GBM_grid_0_AutoML_20181028_212737_mo~ GBM           auc   0.800
##  6 6       GBM_grid_0_AutoML_20181028_212737_mo~ GBM           auc   0.797
##  7 7       GBM_grid_0_AutoML_20181028_212737_mo~ GBM           auc   0.795
##  8 8       GBM_grid_0_AutoML_20181028_212737_mo~ GBM           auc   0.783
##  9 9       GBM_grid_0_AutoML_20181028_212737_mo~ GBM           auc   0.783
## 10 10      GBM_grid_0_AutoML_20181028_212737_mo~ GBM           auc   0.781
## # ... with 40 more rows

• factor_rev() reverses the factor levels.
  • Without fct_rev:

• With fct_rev:

• scales in the ggplot: scales can be fixed (“fixed”, the default), free (“free”), or free in one dimension (“free_x”, “free_y”)

The course example used a deep learning model where only auc and logloss were available metrics. Other models may return those plus mean_per_class_error, rmse, and mse.

# Visualize the H2O leaderboard to help with model selection
plot_h2o_leaderboard <- function(h2o_leaderboard, order_by = c("auc", "logloss"), 
                                 n_max = 20, size = 4, include_lbl = TRUE) {
     # Setup inputs
     order_by <- tolower(order_by[[1]])
     leaderboard_tbl <- h2o_leaderboard %>% as.tibble() %>%
          mutate(model_type = str_split(model_id, "_", simplify = T) %>% .[,1]) %>%
          rownames_to_column(var = "rowname") %>%
          mutate(model_id = paste0(rowname, ". ", as.character(model_id)) %>% as.factor())
     
     # Transformation
     if (order_by == "auc") {
          data_transformed_tbl <- leaderboard_tbl %>% slice(1:n_max) %>%
               mutate(model_id = as_factor(model_id) %>% reorder(auc),
                      model_type = as.factor(model_type)) %>%
               gather(key = key, value = value, 
                      -c(model_id, model_type, rowname), factor_key = T)
          
     } else if (order_by == "logloss") {
          data_transformed_tbl <- leaderboard_tbl %>% slice(1:n_max) %>%
               mutate(model_id = as_factor(model_id) %>% reorder(logloss) %>% fct_rev(),
                    model_type = as.factor(model_type)) %>%
               gather(key = key, value = value, -c(model_id, model_type, rowname), factor_key = T)
          
     } else {stop(paste0("order_by = '", order_by, "' is not a permitted option."))}
     
     # Visualization
     g <- data_transformed_tbl %>%
          ggplot(aes(value, model_id, color = model_type)) + geom_point(size = size) +
          facet_wrap(~ key, scales = "free_x") + theme_tq() + scale_color_tq() +
          labs(title = "Leaderboard Metrics", subtitle = paste0("Ordered by: ", toupper(order_by)),
               y = "Model Position, Model ID", x = "")
     
     if (include_lbl) g <- g + geom_label(aes(label = round(value, 2), hjust = "inward"))
     return(g)}

# Test Model Visualization
automl_models_h2o@leaderboard %>% plot_h2o_leaderboard(n_max = 10, order_by = "auc")

#ggsave("./images/plotH2oLeaderboard.jpg")

The number of facets changes with the model type. Deep Learning models might provide AUC and Logloss. Other models may return rmse, mse, etc.

9.1.1 Gini Coefficient:

\(AUC = (GiniCoeff + 1)/2)\) | \(Gini=(1-AUC)*2\)

The Gini Coefficient is a popular metric on Kaggle, especially for imbalanced class values.

The Gini coefficient is one of many to measure divergence, i.e. the difference in the distribution of two variables in a data set. Values range from -1 to 1 to indicate total divergence one direction or the other, and 0 percent representing no divergence.

h2o.giniCoef(perf_h2o_GLM)# again gives result based on the test data

## [1] 0.74

h2o.giniCoef(myModel_GLM, train = T, valid = T, xval = T)

## train valid  xval 
##  0.80  0.62  0.69

9.1.2 Logloss

Logloss is a metric that measures the class probability from the model against the actual value in binary format (0/1). See this for more information.

Log Loss quantifies the accuracy of a classifier by penalizing false classifications. Minimizing the Log Loss is basically equivalent to maximizing the accuracy of the classifier. A perfect classifier would have a Log Loss of precisely zero. Less ideal classifiers have progressively larger values of Log Loss.

Log Loss heavily penalizes classifiers that are confident about an incorrect classification. For example, if for a particular observation, the classifier assigns a very small probability to the correct class then the corresponding contribution to the Log Loss will be very large indeed. Naturally this is going to have a significant impact on the overall Log Loss for the classifier. The bottom line is that it’s better to be somewhat wrong than emphatically wrong. Of course it’s always better to be completely right, but that is seldom achievable in practice! There are at least two approaches to dealing with poor classifications:

1. Examine the problematic observations relative to the full data set. Are they simply outliers? In this case, remove them from the data and re-train the classifier.
2. Consider smoothing the predicted probabilities using, for example, Laplace Smoothing. This will result in a less “certain” classifier and might improve the overall Log Loss.

h2o.logloss(perf_h2o_GLM)

## [1] 0.3

h2o.logloss(myModel_GLM, train = T, valid = T, xval = T)

## train valid  xval 
##  0.27  0.30  0.32

9.1.3 Confusion Matrix

h2o.confusionMatrix(perf_h2o_GLM)# test data

## Confusion Matrix (vertical: actual; across: predicted)  for max f1 @ threshold = 0.361470825344551:
##         No Yes    Error     Rate
## No     171  13 0.070652  =13/184
## Yes     13  23 0.361111   =13/36
## Totals 184  36 0.118182  =26/220

h2o.confusionMatrix(myModel_GLM)#training data

## Confusion Matrix (vertical: actual; across: predicted)  for max f1 @ threshold = 0.348527593514236:
##         No Yes    Error       Rate
## No     838  51 0.057368    =51/889
## Yes     58 118 0.329545    =58/176
## Totals 896 169 0.102347  =109/1065

Must understand threshold, precision and recall:

• Top Row - Predictions; Left Column - Actual
• Threshold is the value that determines which class probability is a 1 or 0
  • Example: if a 1 is assigned if the probability is above the threshold value - which is a probability
• F1 is the optimal balance between precision and recall. The threshold maximizes F1. However, this is not always the best case. An expected value optimization is required when the costs of false positives and false negatives are known.
  • F1 = ( * (precision * recall))/ (precision + recall) - hence balances precision & recall
• PRECISION - Measures false positives (eg, predicted to leave but actually stay) Precision = TP/(TP + FP) : It detects how frequently the algo over-picks the YES class
• RECALL - Measures false negatives (eg, predicted to stay but actually leave)
  • Recall = TP/(TP + FN) - It provides a metric for under-picking YES
• Because there are different costs associated with FN and FP; FN often cost a company more.
  This is where Expected Value (EV) comes into play - this discussion is coming.

9.4.1 ROC Plot

The ROC Plot compares the TPR (y-axis) against the FPR (x-axis). TPR is the rate which people are leaving correctly identified as leaving. FPR is the rate people that stay are incorrectly identified as leaving.

# Data prep for plotting; # recall numeric to factor must first convert to character
model_metrics_tbl %>% mutate(
     path = str_split(path, pattern = "/", simplify = T)[, 3] %>% as_factor(),
     auc = auc %>% round(3) %>% as.character() %>% as_factor()) %>% 
     
     ggplot(aes(x = fpr, y = tpr, color = path, linetype = auc)) + geom_line(size = 1) + theme_tq() + scale_color_tq()+
     theme(legend.direction = "vertical") + labs(title ="ROC Plot",   subtitle = "Performance of Top Performing Models") +
     annotate(geom = "text", x = 0.6, y = 0.25, label = "TPR - Rate people correctly identified as leaving") +
     annotate(geom = "text", x = 0.6, y = 0.20, label = "FPR - Rate people that stay incorrectly identified as leaving")

9.4.2 Precision vs Recall Plot

Turn precision / recall trade off into performance. This is another multiple model evaluation.

load_model_performance_metrics <- function(path, test_tbl){
     
     model_h2o <- h2o.loadModel(path)
     perf_h2o <- h2o.performance(model_h2o, newdata = as.h2o(test_tbl))
     
     perf_h2o %>% h2o.metric() %>% as.tibble() %>% mutate(auc = h2o.auc(perf_h2o)) %>% 
          select(tpr, fpr, auc, precision, recall)}

# Data prep for plotting; # recall numeric to factor must first convert to character
model_metrics_tbl %>% mutate(
     path = str_split(path, pattern = "/", simplify = T)[, 3] %>% as_factor(),
     auc = auc %>% round(3) %>% as.character() %>% as_factor()) %>% 
     
     ggplot(aes(x = recall, y = precision, color = path, linetype = auc)) + geom_line(size = 1) + theme_tq() + scale_color_tq()+
     theme(legend.direction = "vertical") + labs(title ="Precision vs Recall Plot",  
                                                 subtitle = "Performance of Top Performing Models") +
     annotate(geom = "text", x = 0.6, y = 0.25, label = "Precision - Predicting employees to leave buy actually stay") +
     annotate(geom = "text", x = 0.6, y = 0.20, label = "Recall - Predicting employee to stay but actually leaves")

> Goal is to reduce the under-prediction i.e., employees predicted to stay but actually leave (Recall)

9.5.1 Gain & Lift Tables

Steps:

2. Split the ranked file into 10 sections (deciles). Grouping into cohorts of most likely groups to least likely groups is at the heart of gain/lift. When we do this, we can immediately show that if a candidate has a high probability of leaving, how likely they are to leave.

3. Number of observations in each decile: group = row_number()

4. Number of actual events in each decile: responses = sum(Attrition == "Yes")

ranked_predictions_tbl %>% mutate(ntile = ntile(Yes, n = 10)) %>% group_by(ntile) %>% 
     summarise(cases = n(), responses = sum(Attrition == "Yes")) %>% arrange(desc(ntile))

## # A tibble: 10 x 3
##    ntile cases responses
##    <int> <int>     <int>
##  1    10    22        17
##  2     9    22         8
##  3     8    22         2
##  4     7    22         3
##  5     6    22         2
##  6     5    22         1
##  7     4    22         3
##  8     3    22         0
##  9     2    22         0
## 10     1    22         0

In ntile = 10, 17 people actually left out of 22 employees! 8 of 22 actually left in the next ntile.

• note the use of ntile also from dplyr. We simply group records into 10 groups. Note the 10th group has the highest values for Yes. Makes the continuous variable Yes into a discrete variable using the new groups called ntile.
• dpplyr::row_numer is useful for calculations to simply make the row numbers of a tibble into a column

ranked_predictions_tbl %>% mutate(ntile = ntile(Yes, n = 10)) %>% group_by(ntile) %>% 
     summarise(cases = n(), responses = sum(Attrition == "Yes")) %>% arrange(desc(ntile)) %>% mutate(group = row_number())

## # A tibble: 10 x 4
##    ntile cases responses group
##    <int> <int>     <int> <int>
##  1    10    22        17     1
##  2     9    22         8     2
##  3     8    22         2     3
##  4     7    22         3     4
##  5     6    22         2     5
##  6     5    22         1     6
##  7     4    22         3     7
##  8     3    22         0     8
##  9     2    22         0     9
## 10     1    22         0    10

Now group - the row numbers - now help. The first group has the highest number of employees that actually left. No longer the the ntile column.

Time to calculate the gain and lift.

5. Cumulative responses: Number of cumulative actual events in each decile: cumulative_responses = cumsum(responses)
6. Percentage responses: pct_responses = responses / sum(responses)
7. Cumulative gain: Cumulatively sum the pct_responses. gain = cumsum(pct_responses)
   • Think of this as what we gained by using the model. If we focused on the first 2 groups, we gain the ability to target ~70% of our quitters using the model. Remember, this is measuring employees that actual left and the model is identifying 70% of them!
   • gain is really the cumulative_capture_rate
8. A calculation is needed before we can measure lift. cumulative_pct_cases = cumsum(cases)/sum(cases). This is a simple calculation. Since we have 10 groups, this calculation will start at 0.100 and end with 1.00 step by 0.100 so: .1, .2, 3. . .
9. Calculate lift. Think lift as a multiplier between what we gained divided by what we expected to gain with no model. For example, if we focused on the first 2 decile groups, we gained the ability to target ~70% of our quitter but we expected to be able to target only 20% (2 deciles (cummulative_pct_cases) = 20%)
   • lift is just another way to show the gain
   • Below, lift = 4.7 simply means you are 4.7X more likely to identify quitters using the model

calculated_gain_lift_tbl <-  ranked_predictions_tbl %>% mutate(ntile = ntile(Yes, n = 10)) %>% group_by(ntile) %>% 
     summarise(cases = n(), responses = sum(Attrition == "Yes")) %>% arrange(desc(ntile)) %>% 
     mutate(group = row_number()) %>% select(group, cases, responses) %>% 
     mutate(cumulative_responses = cumsum(responses), pct_responses = responses / sum(responses),
          gain = cumsum(pct_responses),
          # we gain ability to target 75% of quitter focused on 1st 3 groups
          cumulative_pct_cases = cumsum(cases)/sum(cases), lift = gain/cumulative_pct_cases,
          gain_baseline = cumulative_pct_cases, lift_baseline = gain_baseline / cumulative_pct_cases)
calculated_gain_lift_tbl %>% glimpse()

## Observations: 10
## Variables: 10
## $ group                <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
## $ cases                <int> 22, 22, 22, 22, 22, 22, 22, 22, 22, 22
## $ responses            <int> 17, 8, 2, 3, 2, 1, 3, 0, 0, 0
## $ cumulative_responses <int> 17, 25, 27, 30, 32, 33, 36, 36, 36, 36
## $ pct_responses        <dbl> 0.472, 0.222, 0.056, 0.083, 0.056, 0.028,...
## $ gain                 <dbl> 0.47, 0.69, 0.75, 0.83, 0.89, 0.92, 1.00,...
## $ cumulative_pct_cases <dbl> 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0...
## $ lift                 <dbl> 4.7, 3.5, 2.5, 2.1, 1.8, 1.5, 1.4, 1.2, 1...
## $ gain_baseline        <dbl> 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0...
## $ lift_baseline        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1

Above, uncertain why gain_baseline and lift_baseline were added. (BTW, `lift_baseline is always = 1.)

9.5.2 h2o:gainsLift()

h2o.gainsLift streamline the manual gain-lift code above. h2o.performance() is the starting point for the H₂O fain and lift function.

gain_lift_tbl  <-  perf_h2o_GLM %>% h2o.gainsLift() %>% as.tibble()

gain_lift_tbl %>% glimpse()

## Observations: 16
## Variables: 13
## $ group                    <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12...
## $ cumulative_data_fraction <dbl> 0.014, 0.023, 0.032, 0.041, 0.050, 0....
## $ lower_threshold          <dbl> 0.79583, 0.69841, 0.65794, 0.63159, 0...
## $ lift                     <dbl> 6.11, 6.11, 6.11, 6.11, 3.06, 3.89, 2...
## $ cumulative_lift          <dbl> 6.1, 6.1, 6.1, 6.1, 5.6, 4.7, 3.9, 3....
## $ response_rate            <dbl> 1.000, 1.000, 1.000, 1.000, 0.500, 0....
## $ score                    <dbl> 0.8419, 0.7452, 0.6622, 0.6449, 0.625...
## $ cumulative_response_rate <dbl> 1.00, 1.00, 1.00, 1.00, 0.91, 0.77, 0...
## $ cumulative_score         <dbl> 0.84, 0.80, 0.76, 0.74, 0.72, 0.64, 0...
## $ capture_rate             <dbl> 0.083, 0.056, 0.056, 0.056, 0.028, 0....
## $ cumulative_capture_rate  <dbl> 0.083, 0.139, 0.194, 0.250, 0.278, 0....
## $ gain                     <dbl> 511, 511, 511, 511, 206, 289, 122, 12...
## $ cumulative_gain          <dbl> 511, 511, 511, 511, 456, 372, 289, 24...

• group = ntiles by 16: 220 leave /16 ntiles = 13.75 employees/group (manual calculation set ntile = 10)

Only a subset of the H₂O data is needed:

gain_lift_tbl  <-  perf_h2o_GLM %>% h2o.gainsLift() %>% select(group, cumulative_data_fraction, cumulative_capture_rate, cumulative_lift) %>%  as.tibble()
gain_lift_tbl

## # A tibble: 16 x 4
##    group cumulative_data_fraction cumulative_capture_rate cumulative_lift
##    <int>                    <dbl>                   <dbl>           <dbl>
##  1     1                   0.0136                  0.0833            6.11
##  2     2                   0.0227                  0.139             6.11
##  3     3                   0.0318                  0.194             6.11
##  4     4                   0.0409                  0.25              6.11
##  5     5                   0.05                    0.278             5.56
##  6     6                   0.1                     0.472             4.72
##  7     7                   0.15                    0.583             3.89
##  8     8                   0.2                     0.694             3.47
##  9     9                   0.3                     0.75              2.5 
## 10    10                   0.4                     0.833             2.08
## 11    11                   0.5                     0.889             1.78
## 12    12                   0.6                     0.917             1.53
## 13    13                   0.7                     1                 1.43
## 14    14                   0.8                     1                 1.25
## 15    15                   0.9                     1                 1.11
## 16    16                   1                       1                 1

This is all that is needed to create out gain and lift plot.

9.5.3 Gain Plot

• Use gather() to stack the data to use ggplot color, group aesthetics or facet_wrap
• for the gain data, cumumalitive_lift is not needed (but it is later for lift)

gain_transformed_tbl <-  gain_lift_tbl %>% 
     select(group, cumulative_data_fraction, cumulative_capture_rate) %>% mutate(baseline = cumulative_data_fraction) %>% 
     rename(gain = cumulative_capture_rate) %>% gather(key = key, value = value, gain, baseline)
gain_transformed_tbl

## # A tibble: 32 x 4
##    group cumulative_data_fraction key    value
##    <int>                    <dbl> <chr>  <dbl>
##  1     1                   0.0136 gain  0.0833
##  2     2                   0.0227 gain  0.139 
##  3     3                   0.0318 gain  0.194 
##  4     4                   0.0409 gain  0.25  
##  5     5                   0.05   gain  0.278 
##  6     6                   0.1    gain  0.472 
##  7     7                   0.15   gain  0.583 
##  8     8                   0.2    gain  0.694 
##  9     9                   0.3    gain  0.75  
## 10    10                   0.4    gain  0.833 
## # ... with 22 more rows

gain_transformed_tbl %>% ggplot(aes(x = cumulative_data_fraction, y = value, color = key)) +
     geom_line(size = 1.5) + theme_tq() + scale_color_tq() +
     labs(title = "Gain Chart", x = "Cumulative Data Fraction", y = "Gain")

9.5.4 Lift Chart

• Do not need cumulative_capture_rate
• lift-Baseline always = 1

lift_transformed_tbl <-  gain_lift_tbl %>% 
     select(group, cumulative_data_fraction, cumulative_lift) %>% mutate(baseline = 1) %>% 
     rename(lift = cumulative_lift) %>% gather(key = key, value = value, lift, baseline)
lift_transformed_tbl

## # A tibble: 32 x 4
##    group cumulative_data_fraction key   value
##    <int>                    <dbl> <chr> <dbl>
##  1     1                   0.0136 lift   6.11
##  2     2                   0.0227 lift   6.11
##  3     3                   0.0318 lift   6.11
##  4     4                   0.0409 lift   6.11
##  5     5                   0.05   lift   5.56
##  6     6                   0.1    lift   4.72
##  7     7                   0.15   lift   3.89
##  8     8                   0.2    lift   3.47
##  9     9                   0.3    lift   2.5 
## 10    10                   0.4    lift   2.08
## # ... with 22 more rows

lift_transformed_tbl %>% ggplot(aes(x = cumulative_data_fraction, y = value, color = key)) +
     geom_line(size = 1.5) + theme_tq() + scale_color_tq() +
     labs(title = "Lift Chart", x = "Cumulative Data Fraction", y = "Lift")

9.5.5 Explain Gain/Lift Plots

NOTE: The plots used for explanation differ from the calculated plots in the RMD slightly.

Below, doing nothing - not using the model - 25% of employees would be targeted to leave. With the model, 75% employees would be captured - 3X more!

Strategically targeting 25% of high probability customers should yield 75% of potential positive responses!

Lift is a multiplier: How many positive responses would you get at random X Lift = How many would you get with the model? 3X shown below.

Example: How does lift apply to attrition? - Providing stock options strategically to high performers that are at high risk. Could strategically focus on those in the high risk category that are working overtime. Stock options would cost overall less when compared to offering stock options randomly. Reduces cost by 33%!

Always explain in terms of more customers, less churn, better quality, reduced lead times, or anything else executives care about.

10.1.1 Investigating Predictions

This is a preparatory step before getting into the weeds with lime. This should be familiar and serves as a review.

automl_leader %>% h2o.predict(newdata = as.h2o(test_tbl)) %>% as.tibble() %>% head()

## # A tibble: 6 x 3
##   predict    No    Yes
##   <fct>   <dbl>  <dbl>
## 1 No      0.964 0.0359
## 2 No      0.891 0.109 
## 3 No      0.967 0.0330
## 4 No      0.953 0.0471
## 5 Yes     0.210 0.790 
## 6 No      0.937 0.0631

Focus on record 5 since it is predicted to turnover with a high probability. Get more information about this person so we know what actually happened and also get the employees ID so we know who we are talking about.

predictions_tbl <- automl_leader %>% h2o.predict(newdata = as.h2o(test_tbl)) %>% as.tibble() %>% 
     bind_cols(test_tbl %>%  select(Attrition, EmployeeNumber))
predictions_tbl %>% head()

## # A tibble: 6 x 5
##   predict    No    Yes Attrition EmployeeNumber
##   <fct>   <dbl>  <dbl> <fct>              <dbl>
## 1 No      0.964 0.0359 No                   228
## 2 No      0.891 0.109  No                  1278
## 3 No      0.967 0.0330 No                  1250
## 4 No      0.953 0.0471 No                  2065
## 5 Yes     0.210 0.790  Yes                 1767
## 6 No      0.937 0.0631 No                  1308

We now know the employee actually did turnover.

Now find the detail about employee 1767 that did turnover as predicted.

test_tbl %>% slice(5) %>% glimpse()

## Observations: 1
## Variables: 32
## $ Age                      <dbl> 43
## $ Attrition                <fct> Yes
## $ BusinessTravel           <fct> Travel_Frequently
## $ DailyRate                <dbl> 807
## $ Department               <fct> Research & Development
## $ DistanceFromHome         <dbl> 17
## $ Education                <fct> Bachelor
## $ EducationField           <fct> Technical Degree
## $ EmployeeNumber           <dbl> 1767
## $ EnvironmentSatisfaction  <fct> High
## $ Gender                   <fct> Male
## $ HourlyRate               <dbl> 38
## $ JobInvolvement           <fct> Medium
## $ JobLevel                 <fct> 1
## $ JobRole                  <fct> Research Scientist
## $ JobSatisfaction          <fct> High
## $ MaritalStatus            <fct> Married
## $ MonthlyIncome            <dbl> 2437
## $ MonthlyRate              <dbl> 15587
## $ NumCompaniesWorked       <dbl> 9
## $ OverTime                 <fct> Yes
## $ PercentSalaryHike        <dbl> 16
## $ PerformanceRating        <fct> Excellent
## $ RelationshipSatisfaction <fct> Very High
## $ StockOptionLevel         <fct> 1
## $ TotalWorkingYears        <dbl> 6
## $ TrainingTimesLastYear    <dbl> 4
## $ WorkLifeBalance          <fct> Better
## $ YearsAtCompany           <dbl> 1
## $ YearsInCurrentRole       <dbl> 0
## $ YearsSinceLastPromotion  <dbl> 0
## $ YearsWithCurrManager     <dbl> 0

It appears some factors support attrition and some do not. If we want to know why he turned over, use lime.

10.1.2 Lime For Single Explanation

LIME is used to determine which features contribute to the prediction (and by how much) for a single (local) observation

• 2 steps:
  1.Build explainer with lime()
  2. Create Explanation with explain()

h2o, keras and caret have been integrated into lime(). For other packages, use special functions model_type() and predict_model. http://uc-r.github.io/lime

1. To build an explainer, use the ML model and feature distributions (bins) for the training data.

• drop Attrition for the H₂O model
• bin_continuous makes it easy to detect what causes a continuous feature to have a high feature weight in the lime explanation
• quantile_bins tells how to distribute observations within the bins. If TRUE, cuts will be evenly distributed within each bin
• N_bins: 4-5 bins is typically sufficient. If N_bins = 4,bin_cuts` will have 5 cuts. The bins are then the observations that fall in between the cuts.

explainer <- train_tbl %>% select(-Attrition) %>% 
     lime(model = automl_leader, bin_continuous  = TRUE, N_bins = 4, quantile_bins   = TRUE)
#explainer

The output is voluminous and omitted to save pages of printing! It includes how continuous variables are cut. This is important information.

explainer provides much data - be careful before printing.

2. Create explanation with explain

LIME 6 steps for explanation:

1. Given an observation, permute it to create replicated feature data with slight value modifications.
2. Compute similarity distance measure between original observation and permuted observations.
3. Apply selected machine learning model to predict outcomes of permuted data.
4. Select the m number pf features to best describe predicted outcomes.
5. Fit a simple model to the permuted data explaining the complex model outcome with m number features from the permuted data by its similarity to the original observation.
6. Use the resulting feature weights to explain local behavior.

• n_labels: binary classification therefore just one label
• n_features: number of important features to return to explain the record prediction
• n_permutations: the more test/evaluations, the better the explanations (typically). Lime will create 5000 localized linear models for evaluation.
• kernel_width: Affects the LIME linear model (R^2) and therefore should be tuned to ensure you get the best explanations - IMPORTANT
• Attrition is removed because the data argument must match the format the model requires to predict. Since _H2O requires x to be without the target variable, it must be removed below.
• Use lime::explain because explain is a common function used in other packages.

explanation <- test_tbl %>% slice(5) %>% select(-Attrition) %>% 
     lime::explain(explainer = explainer, n_labels = 1, n_features = 8, n_permutations = 5000, 
                   kernel_width = 0.5)
explanation

##       model_type case label label_prob model_r2 model_intercept
## 1 classification    1   Yes       0.79    0.305          0.0699
## 2 classification    1   Yes       0.79    0.305          0.0699
## 3 classification    1   Yes       0.79    0.305          0.0699
## 4 classification    1   Yes       0.79    0.305          0.0699
## 5 classification    1   Yes       0.79    0.305          0.0699
## 6 classification    1   Yes       0.79    0.305          0.0699
## 7 classification    1   Yes       0.79    0.305          0.0699
## 8 classification    1   Yes       0.79    0.305          0.0699
##   model_prediction            feature feature_value feature_weight
## 1            0.492           OverTime             2         0.1783
## 2            0.492 NumCompaniesWorked             9         0.0887
## 3            0.492           JobLevel             1         0.0896
## 4            0.492      MonthlyIncome          2437         0.0799
## 5            0.492   StockOptionLevel             2        -0.0769
## 6            0.492   DistanceFromHome            17         0.0640
## 7            0.492     BusinessTravel             3         0.0569
## 8            0.492     YearsAtCompany             1        -0.0583
##                         feature_desc
## 1                     OverTime = Yes
## 2             4 < NumCompaniesWorked
## 3                       JobLevel = 1
## 4              MonthlyIncome <= 2935
## 5               StockOptionLevel = 1
## 6            13.8 < DistanceFromHome
## 7 BusinessTravel = Travel_Frequently
## 8                YearsAtCompany <= 3
##                                                                                                          data
## 1 43, 3, 807, 2, 17, 3, 6, 1767, 3, 2, 38, 2, 1, 7, 3, 2, 2437, 15587, 9, 2, 16, 3, 4, 2, 6, 4, 3, 1, 0, 0, 0
## 2 43, 3, 807, 2, 17, 3, 6, 1767, 3, 2, 38, 2, 1, 7, 3, 2, 2437, 15587, 9, 2, 16, 3, 4, 2, 6, 4, 3, 1, 0, 0, 0
## 3 43, 3, 807, 2, 17, 3, 6, 1767, 3, 2, 38, 2, 1, 7, 3, 2, 2437, 15587, 9, 2, 16, 3, 4, 2, 6, 4, 3, 1, 0, 0, 0
## 4 43, 3, 807, 2, 17, 3, 6, 1767, 3, 2, 38, 2, 1, 7, 3, 2, 2437, 15587, 9, 2, 16, 3, 4, 2, 6, 4, 3, 1, 0, 0, 0
## 5 43, 3, 807, 2, 17, 3, 6, 1767, 3, 2, 38, 2, 1, 7, 3, 2, 2437, 15587, 9, 2, 16, 3, 4, 2, 6, 4, 3, 1, 0, 0, 0
## 6 43, 3, 807, 2, 17, 3, 6, 1767, 3, 2, 38, 2, 1, 7, 3, 2, 2437, 15587, 9, 2, 16, 3, 4, 2, 6, 4, 3, 1, 0, 0, 0
## 7 43, 3, 807, 2, 17, 3, 6, 1767, 3, 2, 38, 2, 1, 7, 3, 2, 2437, 15587, 9, 2, 16, 3, 4, 2, 6, 4, 3, 1, 0, 0, 0
## 8 43, 3, 807, 2, 17, 3, 6, 1767, 3, 2, 38, 2, 1, 7, 3, 2, 2437, 15587, 9, 2, 16, 3, 4, 2, 6, 4, 3, 1, 0, 0, 0
##   prediction
## 1 0.21, 0.79
## 2 0.21, 0.79
## 3 0.21, 0.79
## 4 0.21, 0.79
## 5 0.21, 0.79
## 6 0.21, 0.79
## 7 0.21, 0.79
## 8 0.21, 0.79

explanation$model_r2[1]

## [1] 0.305

Want R^2 to be as high as possible - optimize by adjusting kernel_width

In this case, kernel_width = 0.5 works relatively well. Trying 0.2 - 1.5 : 0.2 performed best

explanation <- test_tbl %>% slice(5) %>% select(-Attrition) %>% 
     lime::explain(explainer = explainer, n_labels = 1, n_features = 8, 
                   n_permutations = 5000, kernel_width = 0.2)
explanation$model_r2[1]

## [1] 0.315

explanation <- test_tbl %>% slice(5) %>% select(-Attrition) %>% 
     lime::explain(explainer = explainer, n_labels = 1, n_features = 8, 
                   n_permutations = 5000, kernel_width = 1.0)
explanation$model_r2[1]

## [1] 0.312

Explore other tuning parameters

Note we selected 8 features in the explainer so there are 8 records below.

explanation %>% as.tibble() %>% select(feature:prediction)

## # A tibble: 8 x 6
##   feature   feature_value feature_weight feature_desc     data  prediction
##   <chr>             <dbl>          <dbl> <chr>            <lis> <list>    
## 1 OverTime              2         0.184  OverTime = Yes   <lis~ <list [2]>
## 2 JobLevel              1         0.0964 JobLevel = 1     <lis~ <list [2]>
## 3 NumCompa~             9         0.0814 4 < NumCompanie~ <lis~ <list [2]>
## 4 StockOpt~             2        -0.0758 StockOptionLeve~ <lis~ <list [2]>
## 5 Business~             3         0.0773 BusinessTravel ~ <lis~ <list [2]>
## 6 Distance~            17         0.0635 13.8 < Distance~ <lis~ <list [2]>
## 7 MonthlyI~          2437         0.0601 MonthlyIncome <~ <lis~ <list [2]>
## 8 YearsAtC~             1        -0.0557 YearsAtCompany ~ <lis~ <list [2]>

• feature_value is the actual value used. Fro example, this record indicates 9 companies the employee has worked. There are 17 miles between home and work.
• feature_weight indicates the importance. +/- indicates support (positively correlated) or contradict.
• Positive weights support the positive binary classification (Yes or attrition in this case)

kernel_width - R^2 values: 0.3: 0.303 0.4:0.294; 1.0:0.309; 1.5:0.294

In this case (values may vary in RMD), kernel_width = 1 produced the highest R^2 value.

10.1.3 Visualizing A Single Explanation: plot_features()

• ncol = 1 since we only have one record to examine - at this time.
• note there are 3 red bars - the same number of negative feature_weights in the table above.
• Note 4 < NumCompaniesWorked is a label from continuous binning performed earlier
• Red contradict features represent those features that encourage an employee to stay (supports the NO classification)
• GREEN outweigh the RED
• Explanation Fit is the R^2 value

plot_features(explanation = explanation, ncol = 1)

10.1.4 Visualizing Multiple Explanations: plot_explanations()

Modify the slice function for 20 observations by modifying slice.

explanation <- test_tbl %>% slice(1:20) %>% select(-Attrition) %>% 
     lime::explain(explainer = explainer, n_labels = 1, n_features = 8,
          n_permutations = 5000, kernel_width = 0.2)

explanation %>% as.tibble()

## # A tibble: 160 x 13
##    model_type case  label label_prob model_r2 model_intercept
##    <chr>      <chr> <chr>      <dbl>    <dbl>           <dbl>
##  1 classific~ 1     No         0.964    0.272           0.655
##  2 classific~ 1     No         0.964    0.272           0.655
##  3 classific~ 1     No         0.964    0.272           0.655
##  4 classific~ 1     No         0.964    0.272           0.655
##  5 classific~ 1     No         0.964    0.272           0.655
##  6 classific~ 1     No         0.964    0.272           0.655
##  7 classific~ 1     No         0.964    0.272           0.655
##  8 classific~ 1     No         0.964    0.272           0.655
##  9 classific~ 2     No         0.891    0.335           0.804
## 10 classific~ 2     No         0.891    0.335           0.804
## # ... with 150 more rows, and 7 more variables: model_prediction <dbl>,
## #   feature <chr>, feature_value <dbl>, feature_weight <dbl>,
## #   feature_desc <chr>, data <list>, prediction <list>

range(explanation$model_r2)

## [1] 0.239 0.335

#plot_features(explanation = explanation, ncol = 4)

Too much information. The plot is unreadable. Use another method - plot_explanations

plot_explanations(explanation = explanation)

- The No panel indicates the features that support or contradict attrition. Even when there is support for attrition, other factors balance the employee’s experience. - 5, 11, 20 are likely to leave since they are in the Yes panel. - Overtime = Yes is dark red on the No column contradicting a NO response - If use more than 20 cases, it becomes challenging to decipher. The table view might be easier in those cases. Can always dplyr::summarize on the case number to group individual employees

10.1.5 Recreating plot_features() & plot_explanations()

Extend the Lime plots with `theme_tq.

Single case:{#plotlimefeaturetq}

10.1.6 Aside: Glue

Expressions enclosed by braces will be evaluated as R code. Long strings are broken by line and concatenated together. Leading white space and blank lines from the first and last lines are automatically trimmed. Example:

name <- "Fred"; age <- 50; anniversary <- as.Date("1991-10-12") 
glue('My name is {name},', 'my age next year is {age + 1},', 'my anniversary is {format(anniversary, "%A, %B %d, %Y")}.')

10.1.7 plot_features_tq Function

library(glue)# Helps with text

plot_features_tq <- function(explanation, ncol) {
    data_transformed <- explanation %>% as.tibble() %>%
        mutate(
            feature_desc = as_factor(feature_desc) %>% 
                fct_reorder(abs(feature_weight), .desc = FALSE),
            key = ifelse(feature_weight > 0, "Supports", "Contradicts") %>% 
            fct_relevel("Supports"),
            case_text    = glue("Case: {case}"),
            label_text   = glue("Label: {label}"),
            prob_text    = glue("Probability: {round(label_prob, 2)}"),
            r2_text      = glue("Explanation Fit: {model_r2 %>% round(2)}")) %>%
        select(feature_desc, feature_weight, key, case_text:r2_text)
    
    
    data_transformed %>% ggplot(aes(feature_desc, feature_weight, fill = key)) +
        geom_bar(stat = "identity") + coord_flip() + theme_tq() + scale_fill_tq() +
        labs(y = "Weight", x = "Feature") + theme(title = element_text(size = 9)) +
        facet_wrap(~ case_text + label_text + prob_text + r2_text,
                   ncol = ncol, scales = "free")}

New function example:

explanation %>% filter(case %in% 1:4) %>% plot_features_tq(ncol=2)

10.1.8 Multiple cases:

• mutate(order = order_1 * 1000 + order_2) %>% A trick to get one order that preserves the hierarchy of feature (like Monthly_Income) and feature_value (like 9419) in the explanation tibble data.
  • rank simply ranks values from lowest to highest

plot_explanations_tq <- function(explanation) {
     
     data_transformed <- explanation %>% as.tibble() %>%
          mutate(case = as_factor(case), order_1 = rank(feature)) %>% 
          group_by(feature) %>% mutate(order_2 = rank(feature_value)) %>%
          ungroup() %>% mutate(order = order_1 * 1000 + order_2) %>% 
          mutate(feature_desc = as.factor(feature_desc) %>% fct_reorder(order, .desc = T)) %>%
          select(case, feature_desc, feature_weight, label)
     
     data_transformed %>% ggplot(aes(case, feature_desc)) +
          geom_tile(aes(fill = feature_weight)) + facet_wrap(~ label) + theme_tq() +
          scale_fill_gradient2(low = palette_light()[[2]], mid = "white", high = palette_light()[[1]]) +
          theme(panel.grid = element_blank(), legend.position = "right",
               axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
          labs(y = "Feature", x = "Case", fill = glue("Feature Weight"))}

Note glue in labs statement above. With the added return, it stacks Feature and Weight label

plot_explanations(explanation)

plot_explanations_tq(explanation)

11.3.1 Expected Value Baseline w OT: Step 1

• Recall assess_attrition contains several functions:
  • assess_attrition
  • calculate_attrition_cost - This is the one that is needed.
  • count_to_pct
  • plot_attrition
• bind_cols adds columns from test.tbl to the three columns provided by automl_leader %>% h2o.predict(newdata = as.h2o(test_tbl)) %>% as.tibble()
  • Employee Number useful to identify and track employees
  • Monthly Income is used later to estimate the cost of an OT policy change
  • Overtime is needed to determine whether or not the employee is working OT

predictions_with_OT_tbl <- automl_leader %>% h2o.predict(newdata = as.h2o(test_tbl)) %>% as.tibble() %>% 
     bind_cols(test_tbl %>% select(EmployeeNumber, MonthlyIncome, OverTime))
predictions_with_OT_tbl %>% head()

## # A tibble: 6 x 6
##   predict    No    Yes EmployeeNumber MonthlyIncome OverTime
##   <fct>   <dbl>  <dbl>          <dbl>         <dbl> <fct>   
## 1 No      0.964 0.0359            228          9419 No      
## 2 No      0.891 0.109            1278         13577 Yes     
## 3 No      0.967 0.0330           1250         17779 No      
## 4 No      0.953 0.0471           2065          5390 No      
## 5 Yes     0.210 0.790            1767          2437 Yes     
## 6 No      0.937 0.0631           1308          2372 No

Recall calculate_attrition_cost:

calculate_attrition_cost

## function (n = 1, salary = 80000, separation_cost = 500, vacancy_cost = 10000, 
##     acquisition_cost = 4900, placement_cost = 3500, net_revenue_per_employee = 250000, 
##     workdays_per_year = 240, workdays_position_open = 40, workdays_onboarding = 60, 
##     onboarding_efficiency = 0.5) 
## {
##     direct_cost <- sum(separation_cost, vacancy_cost, acquisition_cost, 
##         placement_cost)
##     productivy_cost <- net_revenue_per_employee/workdays_per_year * 
##         (workdays_position_open + workdays_onboarding * onboarding_efficiency)
##     salary_benefit_reduction <- salary/workdays_per_year * workdays_position_open
##     cost_per_employee <- direct_cost + productivy_cost - salary_benefit_reduction
##     total_cost <- n * cost_per_employee
##     return(total_cost)
## }

Some changes we can make:

• annualize the monthly income to calculate the real salary per employee
• capturing net_revenue_per_employee so it can be manipulated
• cost_of_policy_change = 0 because there is no change in the OT policy - keeping everything the same for the baseline
• There is only attrition cost if the employee leaves. The Yes column is the likelihood of leaving (when combined can get an expected attrition cost - coming soon)
• The simplified equation for expected cost:
  • \(Expected Cost = Yes * (AttritionCost + CostOfPolicyChange) + No * (CostOfPolicyChange)\) - Remember Yes and No are probabilities - Because cost_of_policy_change is 0, \(ExpectedCost = Yes * AttritionCost\)

ev_with_OT_tbl <-  predictions_with_OT_tbl %>% 
     mutate(attrition_cost = calculate_attrition_cost(
               n = 1, salary = MonthlyIncome * 12, net_revenue_per_employee = 250000)) %>% 
     mutate(cost_of_policy_change = 0) %>% 
     mutate(expected_attrition_cost = 
                 Yes * (attrition_cost + cost_of_policy_change) + No * (cost_of_policy_change))
ev_with_OT_tbl

## # A tibble: 220 x 9
##    predict    No    Yes EmployeeNumber MonthlyIncome OverTime
##    <fct>   <dbl>  <dbl>          <dbl>         <dbl> <fct>   
##  1 No      0.964 0.0359            228          9419 No      
##  2 No      0.891 0.109            1278         13577 Yes     
##  3 No      0.967 0.0330           1250         17779 No      
##  4 No      0.953 0.0471           2065          5390 No      
##  5 Yes     0.210 0.790            1767          2437 Yes     
##  6 No      0.937 0.0631           1308          2372 No      
##  7 No      0.956 0.0444             18          2661 No      
##  8 No      0.960 0.0403            460          6347 No      
##  9 No      0.963 0.0370           1369          5363 No      
## 10 No      0.959 0.0414           1040          6347 No      
## # ... with 210 more rows, and 3 more variables: attrition_cost <dbl>,
## #   cost_of_policy_change <dbl>, expected_attrition_cost <dbl>

Now calculate the total attrition cost with OT - the baseline case. Simple to to with summarize.

total_ev_with_OT_tbl <- ev_with_OT_tbl %>% summarize(total_expected_attrition_cost_0 = sum(expected_attrition_cost))
total_ev_with_OT_tbl

## # A tibble: 1 x 1
##   total_expected_attrition_cost_0
##                             <dbl>
## 1                        3078803.

11.3.2 Expected Value No OT: Step 2

• Modify test_tbl to change all OT Yes to No. fct_recode makes this easy!

test_without_OT_tbl <- test_tbl %>% mutate(OverTime = fct_recode(OverTime, "No" = "Yes"))
test_without_OT_tbl$OverTime

##   [1] No No No No No No No No No No No No No No No No No No No No No No No
##  [24] No No No No No No No No No No No No No No No No No No No No No No No
##  [47] No No No No No No No No No No No No No No No No No No No No No No No
##  [70] No No No No No No No No No No No No No No No No No No No No No No No
##  [93] No No No No No No No No No No No No No No No No No No No No No No No
## [116] No No No No No No No No No No No No No No No No No No No No No No No
## [139] No No No No No No No No No No No No No No No No No No No No No No No
## [162] No No No No No No No No No No No No No No No No No No No No No No No
## [185] No No No No No No No No No No No No No No No No No No No No No No No
## [208] No No No No No No No No No No No No No
## Levels: No

• Use test_without_OT_tbl to make new predictions using the same code but using a different data file.
• Use rename OverTime so the OT values are side-by-side

predictions_without_OT_tbl <- automl_leader %>% h2o.predict(newdata = as.h2o(test_without_OT_tbl)) %>% 
     as.tibble() %>% 
     bind_cols(
          test_tbl %>% select(EmployeeNumber, MonthlyIncome, OverTime), 
          test_without_OT_tbl %>% select(OverTime)) %>% 
     rename(OverTime_0 = OverTime, OverTime_1 = OverTime1)
predictions_without_OT_tbl %>% head()

## # A tibble: 6 x 7
##   predict    No    Yes EmployeeNumber MonthlyIncome OverTime_0 OverTime_1
##   <fct>   <dbl>  <dbl>          <dbl>         <dbl> <fct>      <fct>     
## 1 No      0.964 0.0359            228          9419 No         No        
## 2 No      0.953 0.0465           1278         13577 Yes        No        
## 3 No      0.967 0.0330           1250         17779 No         No        
## 4 No      0.953 0.0471           2065          5390 No         No        
## 5 No      0.793 0.207            1767          2437 Yes        No        
## 6 No      0.937 0.0631           1308          2372 No         No

Note all the values for Overtime_1 = No. The new overtime values for Overtime_1 will have an impact on the expected value calculation.

11.3.3 Calculate expected value NO OT

• The cost of a policy change will NOT be 0 this time. Instances when OT changes from Yes to No, there is a cost associated with this:
  • \(cost_of_policy_change = OverTime_0 =="Yes" & OverTime_1 =="No" ~ avg_overtime_pct * attrition_cost\)
  • This cost arises primarily from the lost productivity as a result of not working & producing during OT
    • Add new variable avg_overtime_pct. This is a business decision. The goal is to update the calculation to reflect a new policy based on realistic parameters the organization supports.

avg_overtime_pct <- 0.10

ev_without_OT_tbl <-  predictions_without_OT_tbl %>% mutate(attrition_cost = calculate_attrition_cost(n = 1, salary = MonthlyIncome * 12,
                                                            net_revenue_per_employee = 250000)) %>% 
     mutate(cost_of_policy_change = case_when(
          OverTime_0 =="Yes" & OverTime_1 =="No" ~ avg_overtime_pct * attrition_cost, TRUE ~ 0) ) %>% 
     mutate(expected_attrition_cost = Yes * (attrition_cost + cost_of_policy_change) + No * (cost_of_policy_change))
ev_without_OT_tbl %>% select(c(MonthlyIncome:expected_attrition_cost))

## # A tibble: 220 x 6
##    MonthlyIncome OverTime_0 OverTime_1 attrition_cost cost_of_policy_~
##            <dbl> <fct>      <fct>               <dbl>            <dbl>
##  1          9419 No         No                 72979.               0 
##  2         13577 Yes        No                 64663.            6466.
##  3         17779 No         No                 56259.               0 
##  4          5390 No         No                 81037.               0 
##  5          2437 Yes        No                 86943.            8694.
##  6          2372 No         No                 87073.               0 
##  7          2661 No         No                 86495.               0 
##  8          6347 No         No                 79123.               0 
##  9          5363 No         No                 81091.               0 
## 10          6347 No         No                 79123.               0 
## # ... with 210 more rows, and 1 more variable:
## #   expected_attrition_cost <dbl>

11.3.4 Calculate total attrition cost

total_ev_without_OT_tbl <-  ev_without_OT_tbl %>% summarise(
     total_expected_attrition_cost_1 = sum(expected_attrition_cost))
total_ev_without_OT_tbl